> 93 - 2ol8 ? 


NASA CONTRACTOR 
REPORT 


un 

CN 

CM 


■ 



s 



(✓> 

«< 


CASE FILE 
COPY 


AN INTEGRAL-EQUATION SOLUTION 
FOR TE RADIATION AND SCATTERING 
FROM CONDUCTING CYLINDERS 


by J. H. Richmond 


Prepared by 

THE OHIO STATE UNIVERSITY 
Columbus, Ohio 43212 
for Langley Research Center 


NATIONAL AERONAUTICS AND SPACE ADMINISTRATION • WASHINGTON, D. C. • APRIL 1973 



1. Report No. 

2. Government Accession No. 

NASA CR-22^5 


4. Title and Subtitle M INTEGRAL-EQUATION SOLUTION FOR TE RADIATION 

AND SCATTERING FROM CONDUCTING CYLINDERS 

7. Author(s) 


J. H. Richmond 


9. Performing Organization Name and Address 


The Ohio State University 


Electros cience Laboratory 


Dept, of Electrical Engineering 


Columbus, Ohio 43212 


12. Sponsoring Agency Name and Address 


National Aeronautics & Space Administration 

Washington, DC 20546 



3. Recipient's Catalog No. 

\ 


5. Report Date 
April 1973. 


6. Performing Organization Code 


8. Performing Organization Report No. 

TR-2902- 


10. Work Unit No. 


11. Contract or Grant No. 
NGL 36-008-138 


13. Type of Report and Period Covered 
Contractor Report 


14. Sponsoring Agency Code 



16. Abstract 

The piecewise- sinusoidal reaction technique is applied to low frequency radiation and scatter- 
. ing from noncircular cylinders with perfect or imperfect conductivity.- This report presents 
the theory, computer programs and numerical results for these two-dimensional problems with the 
TE polarization. 


17. Key Words (Suggested by Author(s)) 

Antenna analysis, integral equation antenna 
analysis 


18. Distribution Statement 

Unclassified-unlimited 


19. Security Class'd, (of this report) 
Unclassified 


20. Security Classif . (of this page) 

21. No. of Pages 

Unclassified 

81 


22. Price 

$ 3.00 


For sale by the National Technical Information Service, Springfield, Virginia 22151 





















Page Intentionally Left Blank 



CONTENTS 


Page 


I. 

INTRODUCTION 

1 

II. 

THE REACTION TECHNIQUE 

2 

III. 

TRANSVERSE-ELECTRIC STRIP MONOPOLES 

6 

IV. 

THE SINUSOIDAL STRIP DIPOLE 

8 

V. 

CYLINDERS WITH PERFECT CONDUCTIVITY 

15 

VI. 

CYLINDERS WITH FINITE CONDUCTIVITY 

18 

VII. 

THE EXCITATION COLUMN 

18 

VIII. 

FAR-FIELD RADIATION AND SCATTERING 

24 

IX. 

NUMERICAL RESULTS 

28 

X. 

SUMMARY 

38 

REFERENCES 


39 

Appendix 

I 

MAIN COMPUTER PROGRAM 

41 

II 

Subroutine SORT 

48 

III 

Subroutine CDANT 

50 

IV 

Subroutine ZMM1 

53 

V 

Subroutine ZMM2 

55 

VI 

Subroutine ZMM3 

62 

VII 

Subroutine HANK 

66 

VIII 

Subroutine CROUT 

68 

IX 

Subroutine VNAS 

70 

X 

Subroutine VWAS 

71 

XI 

Subroutine VMLS 

73 

XII 

Subroutine CMLS 

74 

XIII 

Subroutine VFF 

76 

XIV 

Subroutine CFF 

78 

XV 

Subroutine CSURF 

79 


i i i 



I. 


INTRODUCTION 


Low-frequency solutions for nonci rcular conducting cylinders 
have been presented in several published papers[l ,2,3, 4, 5, 6, 7,8]. 
Point-matching procedures were employed by Mei and Van Bladel[l], 
Andreasen[2] and Wallenberg and Harrington[5]. In the point-matching 
technique, the contour of the cylinder is divided into segments and 
the integral equation is enforced at the center of each segment. 

Mei and Van Bladel used rectangular pulses as basis functions 
for expanding the surface-current density. They employed trapezoidal- 
rule integration, sampling the integrand only at the center of each 
segment. Andreasen used Simpson 1 s- rule integration and a piecewise- 
parabolic expansion for the current. Wallenberg and Harrington used 
a second-order polynomial for the expansion functions, with dis- 
continuities in the current density at the endpoints of each segment. 

In this report, Rumsey's reaction concept[9] is employed to 
formulate an integral -equation solution for radiation and scattering 
from cylinders with perfect or imperfect conductivity. A piecewise- 
sinusoidal expansion is employed for the current distribution on the 
conducting surface. The complex coefficients in this expansion 
represent samples of the current function. The unknown current 
distribution is forced to have the correct reactions with sinusoidal 
electric test sources located on the conducting surface. (Since the 
test functions are the same as the expansion functions, this is an 
application of Galerkin's method[10].) This procedure generates a 
system of simultaneous linear equations. Numerical solution of this 
system yields a stationary result for the samples of the current 
distribution. Finally the admittance, gain, far-field pattern and 
echo width are determined from the current distribution. 

The new solution appears to be advantageous with respect to 
computational efficiency, convergence, accuracy and generality. The 
current distribution is represented by a continuous function, the 
impedance matrix is symmetric, and the solution satisfies the 
reciprocity and forward-scattering theorems. 

This report considers two-dimensional electromagnetic problems 
involving infinitely-long conducting cylinders. We are concerned with 
the time-harmonic TE case where the field has no z-dependence, E z 
vanishes everywhere, and the time dependence ej wt is understood and 
suppressed. (The z axis is parallel with the axis of the cylinder.) 
With no significant loss of generality, we restrict our attention to 
polygon cylinders. The surrounding medium is free space. The source 
may be an incident plane wave, a parallel magnetic line-source near 
the cylinder, or an axial-slot aperture on the surface of the cylinder. 



With perfect conductivity, the computer programs will handle open as 
well as closed cylinders, arrays of cylinders, and interior as well as 
exterior sources. With finite conductivity, however, the programs are 
restricted to closed cylinders. 

The remaining text presents the detailed theory and some numerical 
results, and the computer programs are listed in the Appendices. 


II. THE REACTION TECHNIQUE 

The reaction concept and its applications have been discussed by 
Rumsey[9], Cohen[ll], Harrington[12] and Richmond[l 3]. 

Consider the exterior scattering problem illustrated in Fig. la. 
(Radiation problems and open surfaces are discussed later.) In the 
presence of a dielectric or conducting body, the impressed electric 
and magnetic currents ( J_-j ,Mj ) generate the electric and magnetic field 
intensities (£,H). For simplicity, let the exterior medium be free 
space . 


From the surface-equivalence theorem of Schelkunoff[14], the 
interior field will vanish (without disturbing the exterior field) 
if we introduce the following surface-current densities 

(1 ) = n x H 

(2) Mg . = E_ x n 

A 

on the closed surface S of the scatterer. (The unit vector n is 
directed outward on S.) In this situation, illustrated in Fig. lb, 
we may replace the scatterer with free space without disturbing the 
field anywhere. 

By definition, the incident field (Ej ,Hj ) is generated by ( Jj ,M-j ) 
in free space, and the scattered field is: 


(3) 


= - " -i 

(4) 

Ms 

-H 


When the surface current (Oc.t^) radiates in free space, it generates 
the field (Es.Hs) in the exterior and (-£j,-H.i) in the interior region. 
This result, illustrated in Fig. lc, is deduced from Fig. lb and the 
superposition theorem. 
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Fig. la. The source (Jj ,Mj ) generates the field (E_,H_) with 
scatterer. 



Fig. lb. The interior field vanishes when the currents (Og.Ms) 
are introduced on the surface of the scatterer. 
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Fig. 1c. The exterior scattered field may be generated by 
(Js,Ms) in free space. 
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With the scatterer replaced by free space, we have noted in 
Fig. lb that the interior region has a null field. As shown in Fig. 
2, we place an electric test source Jj- in this region and find from 
the reciprocity theorem that 


( V s 

\ \ TEST 1 \ 

\ \SOURCE ^ 

.1 


(Js .Ms) ^ 

V's, FREE 
V 


space/ 



FREE SPACE 


Fig. 2. An electric test source is positioned in 
the interior of the scattering region. 


(5) 


U 

<:> 

S 




- W 


ds + 


(ij-it-Ml-Ht) dv=0 


where (it >Ht ) 1S the free-space field of the test source. In words, 

Eq. (5) states that the interior test source has zero reaction with 
the other sources. This "zero-reaction theorem" was developed by 
Rumsey[9], 

Equation (5) is the integral equation for the scattering 
problem, and our objective is to use this equation to determine the 
surface-current distributions and To accomplish this, we expand 
these functions in finite series, so there will be a finite number N 
of unknown expansion constants. Next we obtain N simultaneous linear 
equations to permit a solution for these constants. One such equation 
is obtained from Eq. (5) each time we set up a new test source. 

The magnetic current Ms vanishes if the scatterer is a perfect 
conductor. We assume a finite conductivity and use the impedance 
boundary condition: 


(6) 


% 


= Z, 


0 

-s 


x n 


where Z s denotes the surface impedance. 
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For two-dimensional problems involving cylindrical scatterers, 
(Jj and Ms are functions only of the position z around the contour C 
of the cylinder. If Ji vanishes, Eqs. (5) and (6) yield 


(7) 


o J_ ' [E - (n x H ) Z ] d£ = M. -H ds 
— s -m -m s — i — m 


where (E^H^) denotes the free-space field of test-source m. 

We represent the electric current distribution as follows: 

(8) - J I n 

where the complex constants I n are samples of the function J S U). 

The vector functions J^U) are known as basis functions, subsectional 
bases, expansion functions or dipole modes. We employ expansion 
functions ^ and test sources with unit current density at their 
termi nals . 


From Eqs. (7) and (8) we obtain the simultaneous linear equations 


(9) 

where 


N 

l I Z = V m with m = 1 , 2, 3, • • • N 
“i n mn m 


(10) 


'mn 


J U) • [E - (n x H ) Z ] d£ = 
— n — m — m s 




m 


(ID 


V 

m 



M- -H ds 
—i -m 


] 




d2. . 


In Eqs. (10) and (11) the integrations extend over the region where 
the integrand is non-zero. For example, region n is that portion of 
the contour C covered by the expansion function Jfl. Region m covers 
the interior test source J^. The reciprocity theorem relates the first 
and second integrals in Eq. (10). In the second integral, E^ is the 
free-space field generated by ^ and the associated magnetic current 
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For computational speed and storage, it will be advantageous 
to have a syrrmetric impedance matrix l m . Furthermore, the test 
sources should be selected to yield a well -conditioned set of simul- 
taneous linear equations. For these reasons and to obtain closed 
forms for some of the integrals in Eqs. (10) and (11), we employ 
test sources of the same size, shape and functional form as the 
expansion functions Jn . Finally, we position the interior test 
sources a small distance 6 from surface S and take the limiting form 
of the integrals as 6 tends to zero. 

In this section we have considered explicitely the exterior 
scattering problem. With a slight change in wording, we could make 
the discussion apply equally well to the interior scattering problem. 

To accomplish this, replace "interior region" with "source-free region" 
and replace "exterior region" with "source region". Thus, the unit 
vector n is directed into the source region (which contains M -j ) , and 
we let the test sources approach surface S from the source-free region. 

The next two sections discuss the electric strip dipoles which 
are employed as test sources and expansion modes. Since each dipole 
is comprised of two strip monopoles, the monopoles are considered 
first. 


III. TRANSVERSE-ELECTRIC STRIP MONOPOLES 

Consider the "strip monopole" illustrated in Fig. 3. This 


A 



Fig. 3. An electric strip monopole and 
the coordinate system. 
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source is an electric surface-current density = x J(x) located on 
the xz plane. The source has width h and infinite length and radiates 
in free space. By superposition, the scattered field of a perfectly 
conducting polygon cylinder may be regarded as the field of an array 
of strip monopoles. For the strip monopole shown in Fig. 3, the 
potentials and fields are: 


h 


(12) 

A = 



1 

0 H (kp) dx‘ 
— o 

l 



h 


(13) 

V = 

1 

4ue 

kn 

" ' 8 

O' H Q (kp) dx' 
h 

(14) 

V 

'•J[H (kp) + H 2 (kp) cos (2^) ) 3 dx 
h 


kn 

“ ' 8 

(15) 

E y 

[ J Hp(kp) sin ( 24 ) ) dx 1 • 
0 

(16) 

H = 

h 

_ 

4 

i 

J x p (kp ) dx’ 

) 

where 




07) 

p = 

\|(x - > 

) 2 + / 

(18) 

k = 

W /y£~ 


(19) 

n - 

^y/e 


(20) 

O' ^ 

= dJ 
dx 1 



and the superscript (2) is understood on the Hankel functions H 0 , H-j 
and H 2 - A useful alternative form for the electric field is 

(21 ) £ = - juA - vV 
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h ~ h 

(22) E = -j^ J H Q (kp) dx' + J p J* H^kp) dx' . 

0 0 

For most current functions J(x), the field integrals must be 
evaluated with infinite-series expansions or numerical integration 
procedures. For the sinusoidal current distribution, however, E x is 
obtained rigorously in simple closed form. Thus, if 

(23) J(x) = x [1-j sin (kh - kx) + I 2 sin(kx)] / sin(kh) 
then 

( 2 4 ) E x = 4 S in (kh ) ^1 H o^ kp l^ cos ( kh ^ " ! 1 H 0 ( kp 2 ^ 

+ I 2 H Q (kp 2 ) cos(kh) - I 2 H o (kp-, )] 

where I-j and Io represent J(x) at x = 0 and x = h, respectively. The 
current distribution in Eq. (23) implies line charges at the edges of 
the strip monopole. Since our model of the polygon cylinder will have 
no line charges, the line-charge field contributions are not included 
in Eq. (24). We obtained Eq. (24) by integrating the field of the 
sinusoidal electric line source. Unfortunately, the integral for Ey 
must be evaluated by numerical methods. 


IV. THE SINUSOIDAL STRIP DIPOLE 

A planar strip dipole is illustrated in Fig. 4a. This dipole 
lies in the xz plane and has infinite length in the z direction. The 
surface-current density is 


(25) 

(26) 


. sin k(x - x-j ) 
x sin k(x 2 - x^ 


„ sin k(x 3 - x) 
x sin k(x 3 - x 2 ) 


for x-j < x < x 2 


for x 2 < x < x 3 


As indicated in Fig. 4b, the current density vanishes at the edges 
X] and X 3 , is continuous across the terminals at X 2 and has a slope 
discontinuity at X 2 * 


8 



y 


J j 



Fig. 4a. A planar strip dipole with edges at 
x-j and x 2 and terminals at x 2 - 



X l x 2 *3 


Fig. 4b. The current-density distribution ^ on 
the sinusoidal strip dipole. 


The sinusoidal strip dipole is a hypothetical source in free 
space. The current distribution on a conducting strip is not sinusoidal. 

Figure 5 illustrates a strip V-dipole . Distance along the dipole 
arms is measured by the coordinates s and t with origin at the terminals 
0. The surface-current density is 

. sin k(s 1 - s) 

(27) A = ’ s — iTFTki^ — ' arra s 

. sin k(t, - t) 

( 2S) 0 = t on- arm t 


9 



TABLE I 

Self Impedance of Center-Fed Strip-Dipole Shown in Figure 5 


4 > - 

h/A = 0.05 

h/A = 0.10 

h/A = 0.15 

h/A = 0.20 

45° 

90° 

135° 

180° 

0.11 -j 14.1 
0.38 -j 20.9 
0.64 -j 24.3 
0.75 -j 25.3 

0.47 -j 13.4 
1.59 -j 19.3 
2.68 -j 22.2 
3.11 -j 23.0 

1.20 -j 12.2 
3.94 -j 17.1 

6.49 -j 19.5 

7.50 -j 20.3 

2.53 -j 10.4 
8.10 -j 14.0 
12.95 -j 16.4 
14.77 -j 17.4 



Fig. 5. Nonplanar strip dipole with edges at 
s-| and t-| and terminals at 0. 

/V A 

where the unit vectors s and t are perpendicular to the z axis. Thus, 
the current density vanishes at the edges s-| and t] and has unit value 
at the terminals 0. The edges of the dipole are parallel with the z 
•axis. If the wedge angle \p is adjusted to 180 degrees, the V-dipole 
in Fig. 5 reduces to the planar dipole in Fig. 4-. 

Having defined the sinusoidal strip dipole, we are now in a 
position to explain its relevance. We shall use the dipole current 
distribution (Eqs. (27) and (28)) as the basis function (J^ in Eq.(8)) 
for expanding the unknown current distribution induced on a conducting 
cylinder. Furthermore, strip dipoles will be employed as test sources 
with the reaction concept to solve the integral equation. 
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By superposition, the field of the strip dipole in Fig. 5 is the 
sum of the field contributions from monopoles s and t. The field 
from each monopole can be calculated from Eqs . (12) through (24) with 
the appropriate coordinate transformations. 


Although the sinusoidal strip dipole is a hypothetical source, 
it is useful to define its self impedance with the induced-emf formu- 
1 ati on : 


(29) 


Z 


S 1 

J(s)-E ds 

* 

0 


( J(t)*Edt 
0 


where c[(s) and J^t) are given by Eqs. (27) and (28) and £ is the - 
free-space field of the strip dipole. The reciprocal of Z yields 
the admittance per unit length of the strip dipole. Table I lists 
the self impedance of a center-fed sinusoidal strip dipole as. a function 
of the angle tp and the segment length s-j = t] = h. 

The mutual impedance between two strip dipoles is defined by 


(30) Z ]2 


where ^(s) and J£(t) are given by Eqs. (27) and (28) and Ej is the 
free-space field of the first dipole. 

Figure 6 illustrates a pair of center-fed planar strip dipoles, 
and Table II lists their mutual impedance Z-j 2 - Here p and <j> specify 
the relative positions of the dipoles, 3 specifies the relative 
orientation, and each dipole has the same segment length h. 

Table III lists the mutual impedance Z-j 2 of the overlapping 
strip dipoles shown in Fig. 7. Each of the three segments has the 
same length h, dipoles 1 and 2 share the center segment, and the 
angles are identical. Finally, Table IV gives the mutual impedance 
of the overlapping strip dipoles shown in Fig. 8 . These dipoles share 
an end segment. 


^(s)*^ ds - 


(t ) ' E_i dt 
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TABLE II 

Mutual Impedance of Center- Fed Planar Strip-Dipoles Shown in Figure 6 

Segment length: h/A = 0.1 

Distance between midpoints: p/A =0.3 
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TABLE III 

Mutual Impedance of Center-Fed Overlapping Dipoles Shown in Figure 7 



h/A = 0.05 

h/A = 0.10 

h/A = 0.15 


60° 

90° 

120° 

150° 

180° 

-0.09 +j 8.89 
0.01 +j 6.99 
0.29 +j 6.67 
0.60 +j 6.65 
0.74 +j 6.67 

-0.33 +j 9.85 
0.14 +j 8.26 
1.31 +j 8.22 
2.49 +j 8.16 
2.97 +j 8.10 

-0.59 +j 11.9 
0.77 +j 10.8 
3.52 +j 10.6 
5.88 +j 9.8 
6.73 +j 9.2 

-0.66 +j 15.9 
2.64 +j 15.2 
7.88 +j 13.9 
11.29 +j 11.0 
12.23 +j 9.6 



X 


Fig. 7. Overlapping strip dipoles share 
the center segment. 
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TABLE IV 

Mutual Impedance of Center-Fed Overlapping Dipoles Shown in Figure 8 


a 

h/X = 0.05 

h/X = 0.10 

h/X = 0.15 

h/X = 0.20 

30° 

60° 

90° 

120° 

150° 

0.60 -j 16.4 
0.28 -j 11.7 
0.00 -j 8.25 
-0.09 -j 5.17 
-0.04 -j 2.21 

2.51 -j 14.7 
1 .18 -j 10.7 
0.04 -j 7.82 
-0.38 -j 5.14 
-0.18 -j 2.32 

6.03 -j 12.8 
2.88 -j 9.43 
0.20 -j 6.93 
-0.86 -j 4.84 
-0.46 -j 2.43 

11.83 -j 11.0 
5.80 -j 7.84 
0.71 -j 5.29 
-1.56 -j 3.99 
-0.94 -j 2.50 



Fig. 8. Dipoles 1 and 2 share one segment 
in a Y configuration. 
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V. CYLINDERS WITH PERFECT CONDUCTIVITY 

Consider a perfectly conducting polygon cylinder with contour 
C which may be open or closed. Let J s (&) denote the surface-current 
density induced on the cylinder. If the cylinder is closed, J s will 
flow entirely on the outside or inside of surface S in accordance with 
the location of the source. If the cylinder is open, currents will 
flow on both sides of the thin conducting surface. In all cases, we 
let J s denote the total current density. 

Figure 9a illustrates a perfectly conducting polygon cylinder 
illuminated by a parallel magnetic line source Mj . Let I] and Ig 
denote the current density J s at the comers of the polygon. The 
current J s vanishes at the edges 0 and O'. Let us define two strip- 
dipole mode currents on the cylinder. Mode 1 extends from point 0 to 
point 2 and has terminals at point 1. Mode 2 extends from 1 to 3 
with terminals at 2. Each mode has a sinusoidal current distribution 
and unit terminal current as in Eqs. (27) and (28). Now we represent 
J S U) as the superposition of the two modal currents with weightings 
I] and I 2 . This gives a piecewise-sinusoidal expansion for J S U) 
with two unknown constants I] and I 2 . (In practice we require a 
minimum of around 16 unknowns to obtain accurate results.) 

In the exact solution, the tangential electric field vanishes 
everywhere on contour C. Thus if we move an electric test probe to 
the conducting surface, as in Fig. 9b, the open-circuit voltage at 
its terminals will read zero. To determine N current samples, we 
make N independent probing tests. The probes may be real (thin-wire 
V-di poles) or hypothetical (electric line sources or strip dipoles). 

Now suppose we adjust the currents I n until all the probes 
read zero. This procedure yields a stationary solution for the 
currents I n and, under favorable conditions, tends to the rigorous 
solution as N increases. 

Let Z m denote the mutual impedance between test-probe m in 
Fig. 9b and mode current n in Fig. 9a. The open -circuit voltage in- 
duced in the probe is the sum of the voltage contributions from J s 
and M-:. This voltage must vanish at each probe, leading again to 
Eqs. (9) and (11). With strip-dipole probes, Z mn is given by Eq. (29) 
and Table I for the diagonal elements and by Eq. (30) and Tables II, 

III and IV for the off-diagonal elements. For perfectly conducting 
cylinders, the impedance matrix is symmetric. 

Our piecewise-sinusoidal expansion for the current density J S U) 
satisfies Kirchoff's current law. This follows from the continuous 
nature of each dipole mode in the expansion. In this respect, our 
mode currents resemble the loop currents (as opposed to branch currents) 
in electric circuits. 
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Fig. 9a. Perfectly conducting polygon cylinder with 
parallel magnetic line source NK .. 




Fig. 9b. Electric test-probes 1 and 2 are moved 
to the conducting surface. 
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Fig. 10. Perfectly conducting cylinder illustrates first, 
second and third-order junctions. 


from one comer. This illustrates first, second and third-order 
junctions. (The order of a junction is defined by the number of strips 
extending outward from the junction.) In Fig. 10 the open edge of the 
fin is a first-order junction. Second-order junctions exist at the 
terminals of modes I], I 2 , I 3 and 16- The upper-right corner of the 
square cylinder is a third-order junction. At a junction of order n, 
it is possible to define at least n distinct modes. As a consequence 
of Kirchoff's current law, however, only (n - 1 ) modes are independent. 

In Fig. 10 the current density I 4 flows into the third-order 
junction from below, I 5 into it from the left, and (I 4 + I 5 ) flows out 
of the junction via the fin. 

As illustrated in Fig. 10, a planar conducting surface is 
sometimes divided into two or more segments. The absolute upper 
limit on segment length is a/ 2. For accurate results, however, we 
must use at least four segments per wavelength. 
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VI. CYLINDERS WITH FINITE CONDUCTIVITY 


Section II presents the reaction formulation for closed cylinders 
with finite conductivity. Although the impedance matrix Z mn is symmetri 
for perfectly conducting cylinders, the symmetry disappears with finite 
conductivity. In the square matrix Z mn , the first subscript m denotes 
the row and the second subscript n denotes the column. In the linear 
equations (Eq. (9)), equation m is obtained by enforcing the zero- 
reaction condition with test source m. In Eq. (10), Z m n is the mutual 
impedance between the test source (or test probe) m and the expansion 
mode n. The test sources and the expansion modes are electric strip 
dipoles. With finite conductivity, however, the electric expansion 
modes have associated magnetic currents given by Eq. (6). 

In Eq. (10), the unit normal vector n is directed into the source 
region. By definition, the source region contains the impressed 
current M-j. Thus n reverses direction if M-j is moved across the 
conducting surface S. For perfectly conducting cylinders, the impedance 
matrix Z^ is independent of the source location. With finite con- 
ductivity, however, it is apparent in Eq. (10) that the interior matrix 
differs numerically from the exterior matrix. 

The impedance boundary condition (Eq. (6)) is not rigorous but 
is a reasonable approximation when Z s is small. 

It may be noted from Eqs. (10) and (16) that the mutual impedance 
Z mn between coplanar non-overlapping strip dipoles is independent 
of the surface impedance Z s . 


VII. THE EXCITATION COLUMN 

The complex voltages V m in Eq. (9) form the "excitation column" 
or "excitation vector" in the matrix equation Z mn I n = V m . These 
voltages are independent of the surface impedance Z s . If the impressed 
current is a magnetic line source M ^ , Eq. (11) reduces to 


(31) 




— m —l 


- 

—i -m 


m 


The two forms in Eq. (31) are related by the reciprocity theorem. Both 
forms require numerical integration over test source m. 

If the line source M-s is located at a great distance from the 
cylinder, the incident field (Ej ,Hj ) may be regarded as a plane wave 
wi th 
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(32) 


E. = - ♦ i n H Q e 


jk(x cos <f>.j + y sin<(). ) 


where ^ is the angular coordinate of the source and H 0 is the incident 
magnetic field intensity at the coordinate origin. Figure 11 illustrates 



Fig. 11. A plane wave (E^ ,H^ ) illuminates an 
electric strip dipole. 


an incident plane wave illuminating a strip dipole with edges at points 
1 and 3 and terminals at point 2. Let the sinusoidal electric current 
Jjn fl° w the direction from 1 to 2 and from 2 to 3. The integration 
in Eq. (31) is readily performed to yield 


(33) 


r J ^1 , 

[e - (cos kh, - j cos(a,-<b. ) sin kh, ) e ] 

y = _ n u ! » _J J 

m o k sin kh^ sinCa-j-^) 

[e - (cos kh 9 - j cos(cu-4>.) sin kh„) e ] 

+ L| £ 

n o k sin kh 2 s i n ( ot 2 ~4>-j ) 
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where h] and h2 are the dipole segment lengths and ijq = k(xj cos^-j + 
yj s i n 4> -f ) . The angle between the positive x axis and the vector 
directed to the terminals from point 1 is denoted oq . Similarly a2 
is the angle of the vector directed to the terminals from point 3. 

To investigate the properties of a narrow axial slot in a 
cylinder, we move a magnetic line source Mj to the conducting surface. 
This situation is illustrated in Fig. 12 with the line source located 



Fig. 12. A narrow axial slot is modeled with a magnetic 
line source at the conducting surface. 


at the terminals of dipole mode I2. The arrows on segments (12 and 113 
indicate the reference directions for J2. The circled dot represents 
the magnetic line current Mj flowing in the positive z direction. The 
aperture voltage is Mj with reference polarity indicated by the positive 
and negative signs in Fig. 12. 
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To calculate the excitation column for this situation, we refer 
again to Eq. ( 31 ) and consider the limiting form of the integrals as 
M-j approaches the conducting cylinder from the right. Only the inte- 
grations over segments h£ and (13 require special attention; the others 
are amenable to numerical integration. We position the line source 
M-j at the intersection of segments h2 and I13. Thus the field Ej is 
orthogonal to the conducting surfaces (13 and 113 except on the circular- 
arc region with angle ip. Since the mode current goes sinusoidally 
to zero near the line source, .V] requires integration only over It) 
with no contribution from I12. Likewise V 3 is obtained by integrating 
over h4 with no contribution from h3. The mode current J2 has unit 
value in the vicinity of the line source, and Eq. ( 31 ) yields in the 
limit (as the radius of arc ip tends to zero) 


( 34 ) 



With aperture voltage M-j and aperture current 13, the aperture 
admittance is 


( 35 ) Y = ^/M^ = G + J B mhos/meter. 

As the aperture width tends to zero, the suscepta'nce B is singular. 
The conductance G, however, is well behaved and the time-average 
power delivered by the line source is 

( 36 ) P = |M.j G watts/meter. 

In Fig. 12 , if the line source approaches the conducting 
surface from the left we find 


( 37 ) 


(2t r - ip) M.j 
2 tt 


but the other excitation voltages are not affected. If the cylinder 
is closed, this change in one excitation voltage converts an antenna 
radiation problem to an interior cavity problem. For a perfectly 
conducting cylinder, the aperture conductance will be positive in the 
radiating case and zero in the cavity case. If ip = ir, moving the line 
source across the conducting wall simply reverses the sign of one 
voltage. 

The slot discussed above is called a "one-sided slot" since it 
illuminates only the region on one side of the thin conducting wall. 
Figure 13 a illustrates a one-sided axial slot and the line-source 
model . 
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Fig. 13a. A one-sided slot is modeled with a 
magnetic line source. 



Fig. 13b. A two-sided slot is modeled with a 
magnetic doublet. 


It is convenient to define a magnetic doublet as an array of two 
closely jspaced parallel line sources with equal magnetic currents 
flowing in opposite directions. If the currents M n - are finite and the 
spacing tends to zero, the field Ej vanishes everywhere except in the 
narrow region between the two line sources. As indicated in Fig. 13b, 
the magnetic doublet provides a model for a "two-sided slot" in a 
conducting cylinder. This narrow axial slot radiates into the regions 
on both sides of the thin conducting wall.' In this situation all the 
excitation voltages vanish except V g = Mi , where mode g has the 
doublet at its terminals. This follows from superposition and Eqs. 
(34) and (37). 



An axial slot with finite width is termed a "wide slot" to 
distinguish it from the previously considered "narrow slot" which has 
infinitesimal width. For simplicity assume a uniform electric field 
E a = V/a across the aperture, where V is the aperture voltage and "a" 
is the aperture width. Let the aperture width be divided into an 
integer number of equal segments with length h. The equivalent 
magnetic surface-current density over the aperture is M-j = V/a. For 
a two-sided wide slot, the excitation voltage vanishes for each mode 
outside the aperture region. From Eq. (11), the excitation voltages 
are 


2(1 - cos kh) V 0 

< 38 > v „ - -ran 51 

where e m is unity for the modes in the aperture and zero for those 
outside the aperture. At each edge of the aperture there will exist 
a mode having one segment in the aperture and another segment out- 
side. For these two modes, e m is one-half. 

For a one-sided wide slot the excitation voltages are just one- 
half those in Eq. (38), plus the contributions obtained by numerical 
integration (Eq. (11)) over the segments outside the aperture. 

The complex power supplied by the aperture is 


(39) 


P = VI* = 




where V and I are the aperture voltage and current and the integration 
extends across the aperture. If the electric field is uniform across 
the aperture, E-j = V/a for the two-sided case. It follows that the 
aperture current I is just the average value of J s across the aperture. 
Thus 

N V I | ! • 

(40) I = | JL-H ! ■ ; , 

where V m is given by Eq. (38). 

For the one-sided aperture, E^ = V/(2a) and the aperture current 
I is one-half the aperture -average of J s . In all cases the aperture 
admittance is Y = I/V. 
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VIII. FAR-FIELD RADIATION AND SCATTERING 

The scattered field of a cylinder is the sum of the free-space 
fields generated by the electric surface currents Jg and the magnetic 
surface currents Mg on the cylindrical surface. To obtain the total 
field we add the free-space field (Ej ,Hj ) of the magnetic sources M j . 

For a magnetic line source M z of infinite length located on the 
z axis, the free-space field is 

(41) E = i jkM z H 1 (k P )/4 

(42) H = - z k M z H o (kp)/(4 n ). 

If the magnetic line source is parallel with the z axis and passes 
through the point (x,y), its free-space field at a distant point 
(p ,<p) is 


z k M /2j e"^ kp e jk ^ x cos<i> + y sin<!>) 

(43) H = 

4n Arkp” 


Consider a magnetic surface-current distribution M z (t) on a 
planar strip extending from (x],y|) to (x 2 ,y 2 ) as in Fig. 14. Distance 
from the edge (x-|,y-|) to any point on the strip is measured by the 
coordinate t. For an arbitrary point on the strip. 


(44) x = x-| + t cos a 

(45) y = y i + t sin a . 

From Ea. (43), the free-space field of this source at a distant 
point (p , 4 > ) is 


(46) 


h M z k /ZJ e 


-jkp Jf| 


4n Akp 


M (t) e J ' kct dt 


where 


(47) 4)-| = k(x-| cos<j> + y-j s i n 4> ) 
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(48) 

(49) 

The 

Eq. (16). 
the field 

(50) 


Fig. 14. A planar strip source extends from 
(x-| >y-| ) to (x 2 ,y 2 ). 


\l > 2 = k(x 2 cos <j) + y 2 sin<j>) 
c = cos (<*-<(>). 


free-space field of an electric surface current is given by 
For an electric current J_ = t Jt(t) on the strip in Fig. 14, 
at a distant point (p,<|>) is 


-jkp 

H J = - 1 k *^3" sin(q-j)) e 



j’l'-j 

e 


h 

J (t) e jkct dt. 


It may be noted that Eq. (50) is similar to Eq. (46). 



If the magnetic surface-current density M z (t) in Eq. (46) arises 
through the finite conductivity of the cylinder, the impedance boundary 
condition in Eq. (6) yields 


(51 ) M z (t) = s Z s J t (t) 
where 

(52) s = (t x n)*z = ± 1 


and the unit normal n is directed into the source region. 

From Eqs. (46), (50) and (51) the distant scattered field from 
one segment of the conducting cylinder (the strip in Fig. 14) is 


(53) 
where 

(54) 


/2j (n sin (a-<j>) + s Z ) e 3 p F 
H s = 5 

4n Akp sin kh 


F = k sin(kh) e 1 J t (t) e 3kct dt. 

0 


If the electric current on this segment has a sinusoidal distribution 
as follows: 


(55) J t (t) = 

Eq. (54) yields 


I.| sin(kh-kt) + 1^ sin(kt) 
sin(kh) 


(56) 


1"! j’f'o 

F = — — - — — [e - (cos kh + j c sin kh) e ] 
sin (a-<}>) 

1 2 ■ J'1’2 

+ 2 [e - (cos kh - j c sin kh) e ]. 

sin^(a-<j>) 


In the end-fire directions where (a-<|>) is zero or pi, Eq. (54) yields 



( 57 ) 


j'fh J^p 

F = [e sin kh - kh e jcI-|/2 

j^o j'Pi 

- [e L sin kh - kh e ] jcI 2 /2. 


The field scattered by the cylinder is obtained by summing the contri 
butions from all the segments of the cylinder. 

In plane-wave scattering problems, one is usually interested in 
the echo width W defined as follows: 

(58) W = lim 2 Tip | H s /H 1 | 2 

p-*» 


where H 1 is the incident magnetic field intensity. 

In antenna and radiation problems we are interested in the 
directive gain: 

(59) G = PU)/P av 

where P(^) is the power density in direction <j> and P av is the average 
power density: 


(60) 



1 _ 

2n 


2 7T 

P(<j>) d<j>. 

0 


The power radiated per unit length of cylinder is 


(61) W r = 2-rrp P av . 


For a perfectly conducting cylinder, an alternative expression is 


(62) W p = |V| 2 G 


where V is the terminal voltage and G is the conductance per unit 
length. From Eqs. (59) and (61), 


(63) 


G = 


2irpn I H. 


W.. 
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IX. NUMERICAL RESULTS 


Figure 15 presents the backscattering echo width of circular 



Fig. 15. Backscattering echo width of circular 
cylinder for TE polarization. 


cylinders with perfect and imperfect conductivity. In the reaction 
calculations the cylinder was divided into N segments where N = 12 + 
20 d/x and d is the diameter. Figure 16 shows similar results for a 
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Fig. 16. Broadside backscatter of square 
cylinder for TE polarization. 


square cylinder with 16 segments. Figure 17 illustrates the back- 
scattering characteristics of a circular-sector cylinder. Our data 
in Fig. 17 show excellent agreement with independent calculations by 
Billingsley and Sinclai r[15] . The echo width of the complete circular 
cylinder is 2.224 A. 

Figure 18 shows the conductance of a narrow axial slot in a 
circular cylinder. The number of segments is N = 16 + 20 d/A for 
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ECHO WIDTH/ wavelength 



Fig. 17. Back'scatter echo width of perfectly conducting circular- 
sector cylinder for TE polarization. (Radius: ka = 5) 
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Fig. 18. Conductance of narrow axial slot in perfectly conducting 
circular cylinder for TE polarization. 


diameters up to one wavelength and N = 56 for the larger cylinders. 
Figure 19 shows similar results for a square cylinder with N = 24 + 

20 w/A where w is the width of the cylinder. 

Two methods are available for determining the conductance of a 
narrow axial slot. The simpler and faster method is to take the real 
part of the electric surface -current density at the slot, and the 
alternative is to integrate the far-field power pattern. The two 
methods generally show good agreement if the cylinder is divided into 
an adequate number of segments. For closed noncircular cylinders, 
however, the first method fails when interior modes resonate. With a 
square cylinder, this occurs when the width w is a multiple of x/2. No 
difficulty is encountered in calculating the echo width or the gain 
of a source on or near the cylinder. The conductance data in Fig. 19 
were obtained by far-field integration. 
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CONDUCTANCE ( MILLI MHOS/ WAVELENGTH ) 



Fig. 19, Conductance of narrow axial slot on perfectly 
conducting square cylinder. 


Figure 20 shows the forward gain and backward gain of a narrow 
axial slot in a circular cylinder. 
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Fig. 20. Directive gain of narrow axial slot on perfectly conducting 
circular cylinder for TE polarization. 


Figure 21 illustrates the conductance of a wide axial slot in 
a circular cylinder with 56 segments, and Fig. 22 shows similar data 
for a square cylinder. These conductance data apply when the 
tangential electric field has a uniform distribution across the aperture. 

Figures 23 and 24 show the gain patterns of a magnetic line source 
near a circular and square cylinder, respectively. 
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Fig. 21. Admittance of wide axial slot in perfectly conducting 
circular cylinder for TE polarization. 
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CONDUCTANCE ( MILLIMHOS PER WAVELENGTH ) 






GAIN 



Fig. 23. Gain of magnetic line source near a perfectly conducting 
ci rcul ar cyl inder. 
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GAIN 



Fig. 24. Gain G(<f>) of magnetic line source near a 
perfectly conducting square cylinder. 
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For perfectly conducting strips, our calculations show excellent 
agreement with the eigenfunction solution for broadside backscatter 
and the far-field patterns of one-sided and two-sided narrow axial 
slots. We also have excellent agreement with the eigenfunction solution 
for the far-field pattern of a narrow axial slot in an elliptic 
cylinder[16]. For a magnetic line source near a square cylinder and 
an open circular-arc cylinder, we have excellent agreement with the 
geometrical theory of diffraction. \ 

\ 

\ 

X. SUMMARY \ , 

This report presents the piecewise-sinusoi dal reaction formulation 
for TE radiation and scattering from noncircular conducting cylinders. 
Numerical results are included to show the backscattering echo width of 
circular and square cylinders with perfect and imperfect conductivity. 
Additional data illustrate the admittance and gain of axial-slot antennas 
on circular and square cylinders and the gain of a magnetic line source 
near a cylinder. The computer program and subroutines are presented 
in Appendices. 
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APPENDIX I 

- MAIN COMPUTER PROGRAM 


The main computer program is listed in Fig. 25 in the Fortran 
language. With finite conductivity, the program handles a closed 
cylinder or an array of closed cylinders. With perfect conductivity, 
the program handles open cylinders, closed cylinders, more complicated 
cylinders as in Fig. 10, and arrays of cylinders. By setting the 
integer LOP, the user may select any one of the four excitation options 
listed below. 

1. If LOP = 1, transfer to statement 110 for two-sided narrow 
axial slot. 

2. If LOP = 2, transfer to statement 120 for two-sided wide 
axial slot. 

3. If LOP = 3, transfer to statement 130 for magnetic line 
source or one-sided narrow axial slot. 

4. If LOP = 4, transfer to statement 140 for plane-wave 
scatteri ng. 

Statement 10 reads the first input data defined as follows: 

IWR Positive or negative integer to write or suppress the 
current distribution 

LOP Integer to select option 1, 2, 3 or 4 

NM Number of segments (monopoles) on the contour of the 
cylinder 

NP Number of points on the contour. 

Figures 26 and 27 illustrate typical input data for a square 
cylinder and a strip, respectively. The square cylinder has four 
segments and four points, and the strip has four segments and 5 points. 
(In practice we would divide the contour into a larger number of seg- 
ments for more accurate results.) For a closed cylinder with finite 
conductivity, the points must be numbered in the counterclockwise 
direction (as in Fig. 26) for exterior excitation, and in the clockwise 
direction for interior excitation. For perfectly conducting cylinders, 
the points and segments may be numbered randomly. 
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Next we must give the computer information to indicate which 
points are joined with conducting segments. This data is read in the 
DO LOOP ending with statement 50. IA(J) and I B ( J ) are the index numbers 
of the endpoints of segment J. The segments must be numbered con- 
/ secutively from 1 through NM. The DO LOOP ending with statement 60 

reads the coordinates XC(I) and YC ( I ) of the points in arbitrary 
units. These points, which lie on the contour of the cylinder and 
define its shape and size, must be numbered consecutively from 1 
\ through NP. From XC ( I ) and YC(I), the computer will determine X(I) ' 

\ and Y(I) which denote 2nx/\ and 2iry/\. Statement 80 reads the following 
data : 

CMM Conductivity in megamhos per meter. A negative value 
indicates perfect conductivity. 

DPH Increment in the far-field angle <j> in degrees. 

FMC Frequency in megahertz. 

SCALE Scale factor for multiplying XC and YC to obtain the 

coordinates in meters. A negative value indicates the 
unit of measure for XC and YC is the wavelength. 

TC Thickness of conductor in same units as XC and YC. TC 
is not needed if the cylinder has perfect conductivity. 

For option 1 (two-sided narrow axial slot), statement 110 reads 
the index number IGN of the trade which has the generator at its 
terminals. The modes are set up in subroutine SORT, and statement 55 
writes a list of these modes. Mode I has terminals at the point 12(1) 
and endpoints at 11(1) and 13(1). The reference direction for the 
modal current is from II toward 13, and the reference polarity for 
the modal voltage is negative at II and positive at 13. 

Table V shows the list of modes for the square cylinder and the 


TABLE V 
OUTPUT DATA 


Mode definitions for 
square cylinder in Fig. 26 


I 

WT 

12(1) 

13(1) 

1 

4 

1 

2 

2 

1 

2 

3 

3 

2 

3 

4 

4 

3 

4 

1 


Mode definitions for 


strip in Fig. 27 


I 

11 (I ) 

12(1) 

13(1) 

1 

1 

2 

3 

2 

2 

3 

4 

3 

3 

4 

5 


Yll = .00459 + j .00341 
Gain = 1 .166 at <f> = 0 
Gain = 0.835 at <{> = 180° 
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Scattering cross-section = 0.00163 
Backscatter echo width = 0.00328 





strip. Four modes are defined on the square cylinder, and the terminals 
of mode I are located at point I. Thus to indicate a two-sided narrow 
axial slot at the lower-right corner of the square, we set IGN = 1. 

Three modes are defined on the strip, and the terminals of mode I are 
located at point I + 1. To indicate a slot at the center of the strip 
we set IGN = 2. 

For more complicated cylinders (such as in Fig. 10), it is 
necessary to make a preliminary computer run to obtain the mode 
definitions. Then we can select IGN from the mode list and proceed 
with the final computer run. It is important to note that IGN is a 
mode number and is not always the same as the index number of the 
point where the slot is located. With option 1, the computer will 
write out the slot admittance Yll and proceed to calculate the power 
gain GAIN with increments DPH in the far-field angle PH. 

With option 2 (two-sided wide axial slot), statement 120 reads 
the index numbers JSA and JSB of the first and last segments in the 
aperture. It is assumed that JSB is greater than or equal to JSA 
and that the segments are numbered consecutively across the aperture. 

The computer writes out the admittance Yll and the gain of the wide 
slot. 


For option 3 (magnetic line source or one-sided narrow axial 
slot), statement 130 reads the angle PSI in degrees and the coordinates 
XCS and YCS of the line source in the same units employed for XC and 
YC. PSI is defined in Fig. 12 and Eq. (34) and is used only when the 
line source is located at one of the points (XC,YC) on the contour of 
the cylinder. The computer writes out the admittance Yll and the gain. 
For a unit magnetic line source, the admittance is defined as the 
magnetic field H z at the line source. The singular term in the 
susceptance is suppressed, and the scattered field H| is calculated 
by integrating over the currents Jg and Mg on the cylinder. 

With option 4 (plane-wave scattering), statement 140 reads BSC 
and PHI. Set BSC positive or negative for backscattering or bi static 
data, respectively. For bistatic scattering, PHI is the angle of 
incidence in degrees. The computer writes the scattering angle PH in 
degrees, the echo width EWL in wavelengths, and (in a backscattering 
situation) the extinction cross-section ECS in wavelengths. 

The radiated or scattered power is calculated by numerical 
integration of the far-field power pattern. This result is employed 
to calculate the radiation conductance GR or the scattering cross- 
section SCS immediately below statement 280. These quantities will 
not be accurate, however, unless the angular increment DPH is 
sufficiently small . For a perfectly conducting cylinder, ECS should 
be equal to SCS and GR should be equal to the real part of Yll. 


43 



In some radiation problems the admittance Yll will be inaccurate 
and, as a result, the power gain GAIN will have incorrect normalization. 
This occurs with an axial slot on a closed noncircular cylinder in the 
vicinity of an interior resonance. Even in this case, GAIN is useful 
as the relative far- field power pattern. (The pattern is correct 
even when the normalization is not.) For perfectly conducting 
cylinders, this situation is indicated when GR differs appreciably 
from the real part of Yll. 

Finally, the computer reads JOB and LOP below statement 280. If 
JOB = 10, transfer is made to statement 10 to undertake a new problem. 

If JOB = 80, transfer goes to statement 80 for further calculations on 
a previously considered cylinder with a new conductivity, frequency, 
scale or thickness. If JOB = 300, transfer goes to statement 300 to 
terminate the calculations. If JOB does not equal one of these three 
numbers, transfer goes to statement 110, 120, 130 or 140 according to 
LOP. This is appropriate when we wish to change only the excitation 
of the cylinder. 

In the DATA statement near the beginning of the program, IDM 
denotes the dimensions of the subscripted quantities and INT controls 
the accuracy of certain integrations performed in the subroutines. 
Although INT = 10 is suitable, one may use a smaller or larger value 
to increase the computational speed or the accuracy, respectively. 
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COMPLEX CQT,ZS,Z11 , Z12 ,Y11 ,Y 12 ,HZM ,HZS,HZT 
COMPLEX C (56,56 ),CJ <56 ) ,HJJ( 5 6 ) ,HMM (56 ) ,V.J( 56) 

DIMENSION I A ( 56) ,IB< 56 ) ,11 (56) , I 2 ( 56) , 13(56 ) , JA( 56), JO! 56) 
DIMENSION MD (56,5 ) ,ND (56) ,X ( 56 ) , Y ( 56 ) , D ( 56 ) , XC( 56 > , YC ( 56 ) , DC ( 56 ) 
DATA IDM, INT/56, 10/ 

DATA PI, TP, ETA/3. 14159,6.28318,376.727/ 

2 FORMAT! J.X,8F15.7 ) 

5 FORMAT (1 HO) 

7 FORMAT ( 7F 10 . 5 ) 

8 FORMAT) IX ,1415 ) 

C0T=1 ,41421 4* ETA *CMPLX (l.,-l . ) 

10 READ (5,8 ) IWR ,LOP, NM ,NP 
WRITE(6,8)IWR,L0P, NIT , N P 
WRITEf 6,5 ) 

DO 50 J = 1 , N M 
RFAD(5,8 ) I A ( J ) ,IR (J) 

50 WRITE(6,8)J»IA(J ) , IB ( J ) 

WRITE! 6,5 ) 

CALL SORT ( IA, IB, II , 12, 13 , J A, .»B,MD , ND »NM,NP , N, I DM, MAX, MI N) 
IFIMIN.LT.l .OR. MAX .GT .5 )G0 TO 300 
DO 55 1 = 1, N 

55 WRITE(6, 8)1,11(1) ,12(1), 13(1) 

WR ITE( 6 ,5 ) 

DO 60 1=1, NP 

READ! 5,7) XC ( I ) , Y C ( I ) 

F 1 = I 

60 WRITE! 6,2 )FI ,XC( I ) ,YC ( I ) 

WRITE! 6,5 ) 

DO 70 JAN=1,NM 
K IM= IA( JAN) 

LIM=IB! JAN) 

70 DC ( J AN ) = SORT ( ( XC ( L IM )-XC (KIM)) **2+! YC ( L I M ) -YC ( K I M ) ) **2 ) 

80 READ! 5»7)C MM , DPH ,FMC , S GALE , T C 
WRITE! 6,2 )CMM,DPH ,FMC ,SCALE,TC 
WRITE (6 ,5 ) 

WAVM=300 ./FMC 
TPL =T P 

IF! SCALE. GT .0. )TPL=TP*SCALE/WAVM 
DO 90 I AN= 1 ,N P 
X! IAN) =TPL*XC ( IAN ) 

90 Y( IAN)=TPL*YC ( IAN) 

DO 95 JAN=1,NM 
95 D( JAN ) =TPL*DC (JAN) 

112=1 

ISYM=0 

ZS=(.0,.0) 

TK=TPL-TC 

I F ( C M M . G T . 0 . ) C A L L C S U R F ( C M M , F M C , T K , Z S ) 

IFfCMM.GT.O. ) ISYM--1 

CALL CDANT (C»D,X,Y ,ZS, I A , IB , I 1 » 12,1 3, ISYM 
B , IDM , I NT , JA , J B ,MD,N , ND , NM ,NP ) 

IF! I S YM . E Q • 10 ) GO TO 300 
GO TO ( 110,120,130,140 ), LOP 
110 READ ( 5 , 8 ) I GN 

CALL VNAS ( IDM, IGN, ISYM , IWR , I 12 ,N , C ,C J , Y1 1 ) 

F GN= IGN 

WRITE(6,2)FGN,Y11 
GO TO 200 

120 READ! 5,8)JSA,JSB 

CALL VWAS ( IA, IB, IDM, ISYM, IWR, 11,12,13,1 12, JSA, JSB, MD,N, ND» NM. • 

2 , C ,C J ,D , V J , Y 1 1 ) 

FSA= JSA 


Fig. 25. The MAIN computer program 
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FS8 = J SB 

WRITE! 6t2 ) FSA ,F SB ,Y11 
GO TO 200 

130 READ(5,7)PSI,XCS,YCS 
XS=TPL*XCS 
YS=T PL *YCS 

CALL VMLS ( I A » IB » IDM » INT » IS YM , IWR , 1 1, 1 2 , I 3, 1 12 » MD, N»ND »NM, 
2C,CJ,D,PSI,VJ,X,Y,XS,YS,Y11,ZS) 

WRITE(6,2)PSI,XCS,YCS,Y11 
GO TO 200 

140 READ! 5,7>BSC,PHI 
WRITE! 6,2 )BSC, PHI 
200 WR ITE ( 6 ,5 ) 

IF(L0P.NE.4)G=REAL (Yll) 

INC=— 1 

IF(L0P.EQ.4)INC=1 

IPA=2 

IF(L0P.E0.4)IPA=1 

IF(L0P.NE.4)BSC=-1. 

NPH=360./DPH+1.5 
GR= .0 

DC 280 IPH=IPA,NPH 
FPH= I PH-2 
PH=DPH*FPH 
I F ( IPH.EO.l ) PH=P HI 

CALL VFF ( IAt IB, INC , IDM, ISYM, IWR, I 1, 12, 13, 1 12, LOP , MD , N, NO, NM, 
2C ,CJ,D,EWL,G,GAIN,HJJ ,HMM,HZS,HZT ,PH,ECS,VJ,X,Y,XS,YS,ZS) 

IF ( L OP . N E .4 ) WR I TE (6,2 )PH,GAIN 
I F( L0P.E0.4 ) WR IT E(6,2 )PH,EWL ,ECS 
INC=— 1 

1F( BSC .GT .0.) INC =1 
280 GR=GR+CABS (HZT )**2 
WRITE(6,5) 

SC S= .0174533 *DPH*GR 

I F ( L0P.E0.4)WRITE(6,2)SCS 

GR=ETA*SCS 

I F ( L0P.NE.4 ) WRITE* 6, 2 )GR 
WR ITE( 6 » 5 ) 

READ! 5,8IJOB,LOP 
IF! JOB.EO.IOIGO TO 10 
IF! J0B.E0.80 )G0 TO 80 
IF! JOB.E0.300)GO TO 300 
GO TO! 110, 120, 130, 140), LOP 
300 CONTINUE 
CALL EXIT 
END 


Fig. 25. (Continued) 
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Card 

Number 


1 

-13 4 4 

2 

1 2 

3 

2 3 

4 

3 4 

5 

4 1 

6 

.5 -.5 

7 

.5 .5 

8 

-.5 .5 

9 

-.5 -.5 

10 

-1. 180. 

11 

90. .5 

12 

300 1 



300. .1 .001 

-.5 


Fig. 26. Input data for square cylinder with one-sided narrow 
axial slot. The diagram shows the numbering system 
for the four points on the cylinder. 


Card 







Number 







1 

-1 

4 

4 5 




2 

1 

2 



y 

i 


3 

2 

3 



i 


4 

3 

4 


• • — 

— • — 

— • -9 

5 

4 

5 


1 2 

3 

4 5 

6 


-.5 

.0 




7 


-.25 

.0 




8 


.0 

.0 




9 


.25 

.0 




10 


.5 

.0 




11 


-1. 

90 . 

300 . 

.1 

.001 

12 


-1. 

90 . 




13 

300 

1 






Fig. 27. Input data for bi static scattering by strip. The 
diagram shows the numbering system for the five 
points on the strip. 
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APPENDIX II 
Subroutine SORT 


Subroutine SORT, listed in Fig. 28, defines a set of strip- 
dipole mode currents on the cylinder. The input data IA, IB, NM, 

NP and IDM have been defined already. The output data are defined 
as follows: 

N 

11 ( 1 ) 

12(1) 

13(1) 

JA(I) 

JB(I) 

ND(J ) 

MD(J,K) 

MAX, MIN 

At completion of the DO LOOP ending with statement 20, NJK denotes 
the number of segments intersecting at point K, and JSP is a list of 
these segments. In the DO LOOP ending with statement 22, the computer 
sets up the appropriate number MOD of dipole modes with terminals at 
point K. The first mode is set up on segments JSP (1 ) and JSP (2 ) , the 
second mode on JSP (1 ) and JSP (3) , etc. Since the segments are listed 
in ascending order in JSP, all the modes with terminals at point K 
will share the lowest-numbered segment JSP ( 1 ) having an endpoint at K. 

LAX denotes the largest number of segments intersecting at any 
point, and LIN denotes the smallest. 


Number of dipole modes on the cylinder 

Endpoint on dipole I 

Terminal point on dipole I 

Endpoint on dipole I 

Segment of dipole I 

Segment of dipole I 

Number of dipoles sharing segment J 

A list of the dipoles sharing segment J 

Extreme values of ND ( J ) . 


SUBROUTINE SORT ( TA , IB , II , I 2 , 1 3 , J A , JB , HD, ND, NM , NP , N , I DM , MAX, M IN ) 
DIMENSION JSP (10) ,MD ( IDM, 5 ) ,MD( IDM) 

DIMENSION I A ( IDM ), IB ( IDM) , II ( IDM ) , 12(1 DM) ,I3( IDM) , JA ( IDM) , JBC IDM) 
1=0 
LAX=0 
L I N= 100 
DO 24 K=1,NP 
NJK=0 

DO 20 J=1 »NM 

IND= ( IA(J)-K)*(IR (J)-K) 

I F ( IND.NE.O )G0 TO 20 
N J K= NJK+ 1 
JSP ( NJK )=J 
20 CONTINUE 

IP( NJK.GT .LAX )LA X=NJK 
IF( NJK.LT .L IN)L IN =N JK 
MOD=N JK— 1 

I F ( MOD .L E «0 ) GO TO 24 

Fig. 28. Subroutine SORT. 
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JA 1 = JSP ( 1 ) 

I 1 1= I A { JAI ) 

IF l 1A< JAIJ.EO.K) II1=IB ( JAI ) 

DO 22 IMD=1 f MOD 
1 = 1 + 1 
NJK= IMD+1 
JA ( I ) = JAI 
J8I=JSP(NJK) 

JB ( I ) = JB I 

nm = iii 

1 2 ( I )=K 
1 3 ( I ) = I A ( JBI ) 

IF { IA( JBI J.EQ.K) 13(1 > = IB(J8I ) 

22 CONTINUE 
24 CONTINUE 

I F ( LAX.NE .2 .OR. LIN.NE.2JGQ TO 26 
111 = 11 ( 1 ) 

II ( 1) — 13(1) 

1 3 ( 1 > = 1 1 1 
JA 1= JA ( 1 ) 

JA ( 1 ) = JB ( 1 ) 

J 8 ( 1 ) = JA 1 
26 N= I 

DO 30 J=1,NM 
ND ( J ) =0 
DO 30 K= 1 , 5 
30 MD( J » K ) =0 
DO 40 1=1, N 
J = J A ( I ) 

DO 38 1=1,2 
ND( J ) =ND ( J ) +1 
K= 1 
M=0 

32 MJK=MD(J,K) 

I F ( MJK.NE.O )G0 TO 34 
M= 1 

MD( J,K)=I 
34 K=K+ 1 

I F ( K .GT .5 )G0 TO 38 
I F ( M . E 0 . 0 ) GO TO 32 
38 J = JB( I ) 

40 CONTINUE 
M IN= 100 
MA X= 0 

DO 46 J=1 ,NM 
NDJ=ND( J ) 

I F ( NDJ.GT.MAX )MAX=NDJ 
46 IFtNDJ.LT .MINIM IN =NDJ 
RETURN 
END 


Fig. 28. (Continued) 



APPENDIX III 
Subroutine CDANT 


Subroutine CDANT, listed in Fig. 29, sets up the impedance 
matrix Zy defined by Eqs. (9) and (10) and denoted in the computer 
program by C(I,J). Most of the input data have been defined already, 
the new items are defined as follows: 

D(J) kdj where dj is the length of segment j and k = 2u/x 

X(I),Y(I) kx-j and ky-j where (x-j ,y-j ) are the coordinates of point i 

ZS surface impedance in ohms 

ISYM zero or one for perfectly conducting or finitely con- 

ducting cylinder, respectively. 

The surface impedance ZS is assumed to be uniform over the entire 
conducting surface. If the cylinder has perfect conductivity, Zy is 
syrmietric and CDANT sets up just the upper-right triangular portion of 
C(I,J). With finite conductivity, the entire square matrix is generated. 
If any segment has a length exceeding 0.477 x, CDANT returns with 
C(I,J) = 0 and ISYM = 10. 

C(I , J ) denotes the mutual impedance between the electric test 
dipole I and expansion mode J. To calculate Zy from Eq. (10), we 
integrate over dipole J using the current Jj of dipole J and the 
field of dipole I. It is convenient to express Z-jj as the sum of an 
integral over segment JA(J) and an integral over segment JB(J). (These 
are the two segments of dipole J.) Furthermore, the field of dipole 
I is the sum of the fields from segments JA(I) and JB(I). Thus, Zy 
can be expressed as the sum of four double integrals: 

< 64 > z n-i n U • } i* H <> - p " 4 ds dt 

L K L K 

where Jj is the current density of dipole I. Equation (64) follows 
from Eqs. (10) and (22) for a perfectly conducting cylinder. The 
argument of H 0 and H-] is kp where p is the distance between a source- 
point on- segment K and a point on segment L, and p is the corresponding 
unit vector. The dipole current distributions J-j(s) and Jj(t) are 
sinusoidal as in Section IV and have unit value at the dipole terminals. 
From Eq. (64), it is convenient to consider the dipole-dipole im- 
pedance Zjj to be the sum of four segment-segment impedances Z|_k. 

CDANT selects a segment K (at statement 30) and another segment 
L. As in Eq. (64), K is a segment of test-dipole I and L is a segment 
of expansion mode J. The segment-segment impedance Z[_k is obtained by 
calling a subroutine. In statement 168, this impedance is then added 
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In CDANT, segment K has endpoints KA and KB and segment L has 
endpoints LA and LB. It is convenient to think of KA and KB as points 
1 and 2 on segment K, and LA and LB as points 1 and 2 on L. Now we 
can define four segment-segment impedances P(IS,JS). The first sub- 
script IS refers to the terminal point on segment K, and the second 
subscript JS refers to the terminal point on L. Thus IS = 1 or 2 
if dipole I has its terminal point 12(1) at KA (point 1) or KB 
(point 2), respectively. Similarly JS = 1 or 2 if mode J has its 
terminal point 12 (J ) at LA or LB, respectively. The impedances 
P ( I S , JS ) are defined with the following reference directions for 
current flow: from point 1 toward point 2 on each segment. If dipole 

I has this same reference direction on segment K we set FI = 1 ; 
otherwise FI = -1. Similarly, FJ - ±1 in accordance with the reference 
direction for mode J on segment L. In statement 168 P ( IS, JS) is 
multiplied by FI and FJ before its contribution is added to Z^j. 

CDANT calls subroutine ZMM1 if the segment numbers K and L are 
identical, ZMM2 if segments K and L intersect, or ZMM3 otherwise. 
Subroutine ZMM2 calculates the impedances Q(KK,LL) which are like the 
P(IS,JS) but have different conventions on reference directions and 
subscript meaning. The transformation from the Q impedances to the 
P impedances is accomplished in statement 98. 


SUBROUTINE CDANT! C , D , X , Y , ZS , I A , I B 1 1 1, I 2 , I 3, I-SYM 
2, IDM t INT,JA,JB,MD,N,NDtNM,NP) 

COMPLEX ZS,Pll,P12,P2l,P22,011,Q12,Q21tQ22,P(2,2>,Q!2»2> 
COMPLEX C( I DM , IOM) 

DIMENSION X { I DM ), Y ( I DM) »D ( I DM ) , I A( I DM ) , IB( I DM) , JA ( I DM) , JB( I DM) 
DIMENSION II ( I DM ) , 12 ( I DM-) , 13 ! I DM ) , MD ( I DM , 5 ) < ND I I DM ) 

DO 20 1=1, N 
DO 20 J=1,N 

20 C( I,J)=( .0,.0) 

DMAX= .0 
DO 25 J=1 ,NM 
DK=D( J) 

25 I F ( DK «GT . DMAX ) DMA X=DK 
IF! DMAX.LT .3. )G0 TO 30 
I S YM= 10 
RETURN 

30 DO 200 K= 1 »NM 
NDK=ND ! K ) 

KA= I A ( K ) 

KB= IB ( K ) 

DK = D ( K ) 

DO 200 L=1,NM 
NDL=ND( L ) 

LA= IA ( L ) 

LB= IB( L ) 

DL=D( L ) 

N IL = 0 


Fig. 29. Subroutine CDANT. 
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00 200 11=1, NDK 

1 = MO ( K , 1 1 ) 

F 1= 1 . 

IF(KB.E0.I2(I) ) GO TO 36 
IFtKB.EO.IKI ) )F I=-l. 

IS = 1 

GO TO 40 

36 IF ( KA .EG • 13 ( I ) ) FI =— 1 • 

I S= 2 

40 DO 200 J J= 1 ,NDL 
J=MD( L , JJ ) 

I F ( ISYM.NE.O)GO TO 42 
I F ( I .GT .J ) GO TO 200 
42 F J=1 . 

IF(LB.E0.I2(J ) )G0 TO 46 
IFILB.EO.Il (J) )FJ =-l. 

JS= 1 

GO TO 50 

46 I F ( LA.E0.I3U ))FJ=-1. 

J S=2 

50 I F ( NIL »ME «0 )G0 TO 168 
N IL= 1 

IFIK.EO.L )GQ TO 120 

IND= ( LA-KA )*(LB-K A )*(LA-KB )*t LB-K3) 

I F ( IND.EO.O )G0 TO 80 
C SEGMENTS K AND L SHARE NO POINTS 

CALL ZMM3(X(KA)»Y(KA),X(KB)»Y(KB),X(LA),Y(LA)»X(LB)»Y(LB)»ZS» 
2DK,DL, INT,P(1,1 ), P ( 1,2) ,P(2, 1) ,P( 2,2) ) 

GO TO 168 

C SEGMENTS K AND L SHARE ONE POINT (THEY INTERSECT) 

80 KG= 3 
JM=KB 
JC= KA 
K F=— 1 

IND=( KB— LA )* (KB— LB ) 

I F ( I NO • N E • 0 ) GO TO 82 

JC= KB 

KF = 1 

JM= KA 

KG=0 

82 LG= 3 
J P = L A 
L F= — 1 

IF ( LB .EO.JC )G0 TO 83 
JP= LB 
L F = 1 
LG= 0 

83 SGN=KF*LF 

CALL ZMM2(X( JM ) , Y( JM ) , X ( JC ) , Y ( JC),X( JP ),Y( JP ) , ZS, OK, DL, 

2INT, 0(1,1 ) ,0 ( 1 , 2 ) ,0(2,1 ),Q(2,2)> 

DO 98 KK=1 ,2 
KP= IABS ( KK-KG) 

DO 98 LL=1 ,2 
L P= IABS ( LL-LG) 

98 P(KP,LP)=SGN*Q(KK,LL) 

GO TO 168 

C K=L (SELF REACTION OF SEGMENT K ) 

120 CALL ZMM1 ( DK ,ZS ,P ( 1 ,1 ),P ( 1,2 ) ) 

P(2,1)=P(1»2) 

P( 2,2)=P(1,1 ) 

168 C( I ,J)=C(I,J)+FI *FJ*P ( IS , JS ) 

200 CONTINUE 
RETURN 
END 

Fig. 29. (Continued) 
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APPENDIX IV 
Subroutine ZMMl 


ZMMl, listed in Fig. 30, calculates the self-impedances of a 
sinusoidal strip monopole as defined in Appendix 3. (The terms 
"monopole" and "segment" are used interchangeably.) These impedances 
are obtained in functional form (as opposed to numerical integration) 
with the aid of Eq. (24) and the following: 

h 

» 

(65) k H Q (kx) sin(kx) dx = 2j/u + kh [ H q ( kh ) sin(kh) - H^(kh) cos(kh)] 

0 

h 

(66) k | H Q (kx) cos(kx) dx = kh [H q ( kh ) cos(kh) + H^(kh) si n (kh ) 3 

0 


where k and h are positive. Since ?21 = P -| 2 and ?22 = Pll > ZMMl 
calculates only Pn and P-^. 

The impedance contribution aP arising from finite conductivity 
is 



where the subscripts m and n refer to the test-source and the expansion 
mode, respectively. To explain the factor of one-half in the last form 
of Eq. (67), consider the free-space field of the test monopole. If 
the observer approaches the surface of the monopole, the field is H™ = 

± £ Jm/2. 
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SUBROUTINE ZMM1 ( OK , ZS , P 1 1 , P 12 ) 

COMPLEX ZS,H0,H1,P11,P12 
DATA PI/3.14159/ 

CDK=COS ( DK ) 

SDK = S IN (DK ) 

CALL HANK (DK,HO,Hl ,2 ) 

SDKS=SDK**2 

CDKS=CDK**2 

Pll =-2.*Hl*CDK + HO*SDK+2 .*(.0,1. )*{ l.+CDKS) /PI/DK 

PI 2 =-HO*CDK*SDK+Hl*( l.+CDKS) -4. *(.0,1. )*CDK/PI/DK 

P11=15.*DK*P11/SDKS 

P12=15.*DK*P12/SDKS 

RS= REAL ( ZS ) 

IF(RS.LE.O. )G0 TO 100 
CST=16.*PI*SDKS 
TDK=2.*DK 
CTDK=COS ( TDK ) 

STDK=SIN(TDK) 

P11=P11+ZS*(TDK— STDK )/CST 

PI 2 = PI 2 +ZS*( (l.-CTDK)*SDK + (STDK-TDK)*CDK)/CST 
100 CONTINUE 
RETURN 
END 


Fig. 30. Subroutine ZMM1 . 
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APPENDIX V 
Subroutine ZMM2 


Subroutine ZMM2, listed in Fig. 31, calculates the segment- 
segment impedances of two intersecting strip monopoles. The input 
data are: 

XI, Y1 Coordinates (kx-j ,ky-| ) of the free end of the source segment 

X2, Y2 Coordinates of the point of intersection 

X3, Y3 Coordinates of the free end of segment 2 

DK1 Length kh-j of source segment 

DK2 Length kh£ of segment 2. 


Figure 32 illustrates the intersecting segments after a co- 
ordinate rotation and translation. Point 1 is at the origin and 
point 2 lies on the positive x axis at x = h-j. The new coordinates of 
point 3 are denoted by XB and YB. Segment 1 lies on the x axis and 
a denotes the angle between segment 2 and the x axis. AL, CAL and 
SAL denote a, cos a and sin a. If a vanishes, the two segments are 
coplanar. In this event, the impedances are calculated in closed form 
just above Statement 20. In the coplanar situation the magnetic field • 
of segment 1 vanishes over segment 2 and there is no coupling via 
the surface impedance. 

Figure 32 shows the reference directions for the current density 
^ as used in the definition of the impedances Qll, Q12, Q21 and Q22. 
Figure 32 also shows the endpoint numbering system for segments 1 and 
2 as used in defining these impedances. As usual, the first and second 
subscripts on Q indicate the terminal points on segments 1 and 2, 
respectively. 


If the cylinder has perfect conductivity, the mutual impedance 
between monopoles 1 and 2 is 


(68) Z = 



J 2 (t)-E_ 1 (t) dt 


where J? is the expansion -mode current density on monopole 2 and E_-j 
is the free-space field of monopole 1. Coordinate t denotes the dis- 
tance from the junction to an arbitrary point on segment 2. 

Let E x and E y denote the components of £] . It is convenient 
to let Z = Z x + Zy'where 
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(69) Z x = - cos (a) | J(t) E x (t) dt 

(70) Z y = - sin (a) | J(t) E y (t) dt 


ZMM2 uses Eq. (24) for E^. Comparison of Figs. 3 and 32 shows that 
P2 = t. Thus, one term in Z x can be integrated in closed form via 
Eqs. (65) and (66). This is accomplished immediately after statement 
20, and the result is denoted by SI 1 , SI 2, etc. The remaining term 
in Z x requires the integration of J(t) H 0 (kp-j ) . This is evaluated 
with Simpson's rule in the DO LOOP ending at statement 90, and the 
result is immediately lumped into SI 1 , SI 2 , etc. 

We still have to calculate Z y and the impedance contribution 
aZ from the surface impedance ZS. ^For Z y we obtain E y from the second 
integral in Eq. (22). For aZ we use the^ first expression in Eq. (10) 
and H_ from Eq. (16). Both Z y and aZ have the following form: 

h 2 h l 

(71) R = sin (a) F ] (s ) J 2 (t) (t/p ) ^ (kp ) ds dt 

where 

(72) p 2 = s 2 + 2st cos (a) + t 2 

and F](s) represents J-](s) or J-j(s). Coordinate s denotes the distance 
from the junction to an arbitrary point on segment 1. In Eq. ; (71) 
the integrand is singular at (s,t) = (0,0). To remove the singularity, 
we make three successive changes of variable as follows: 

(73) p = (t + s)//T q = (t - s)//2" 

(74) u = vT p cos (a/2) v = SE q sin (a/2) 

(75) u = p cos tj) v = p sin<j>. 

In the final transformation (Eq . (75)), the integration variable <|> 
should not be confused with the angle <j> shown in Fig. 3. From Eqs. 
(71)- (75), 
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(76) s - p s 1n <“/? - -. ti . 
v ' H Sina 


t = sin (a/2 + 4>) 


Sina 


a/2 p m 


(77) R = 


F.| (s) J 2 ( t ) t ( k P) d P d< f> 


-a/2 0 


h-| sin(a) 

(78) - s -j n ( a /2 - (j,)’ 




h^ sin(a) 

p m _ sin (a/2 + <j>) 5 


* 1*0 


(h ? - h 1 ) tan (a/2) 

t 80 ) tan h-nr, 


ZMM2 evaluates Eq. (77) with Simpson's rule. The p integration is 
performed in the DO LOOP ending at statement 100, and the <j> integration 
is in the outer DO LOOP ending at statement 200. 
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SUBROUTINE ZMM2fXl ,Y1 , X2 , Y2 ,X3 , Y3 ,Z S 
2, DK1, OK 2, I NT, Oil ,012,021,022 ) 

COMPLEX HO , HI , HHO ,HH1 ,SH0,SH1 ,0 1 1 ,0 12 ,Q21 ,022 
COMPLEX OHHO ,DHH 1 , DHO ♦ DH1 , DSHO , OSH 1 

COMPLEX S11,S12,S21,S22,T11,T12,T21, T22 , Y 1 1 , Y 12, Y21 ,Y22 
COMPLEX Dili ,DT12,DT21 ,DT2 2 , DY l 1 ,D Yl 2, DY2 1 » DY22 
COMPLEX ZS ,RKH1 ,SX1 ,SX2,CCP,FUN,CQT 
DATA CCP/ ( .0, .63662 )/ 

DATA PI/3.14159/ 

SDK 1= S I N ( OKI ) 

S0K2=S IN ( DK2 ) 

CDK 1=C0S (OKI) 

CDK2=C0S (DK2) 

CBET=< X2-X1 J/DKl 

SB ET= ( Y2-Y1 )/DKl 

XB= ( X3-X1 )*C8ET+ (Y3-Y1 )*SBET 

YB=-( X3-X1 ) *SB ET+ ( Y3-Y1 )*CBET 

CAL= ( XB-DK1 ) / DK2 

SAL=ABS ( YB/DK2 ) 

CALL HANK (DK2 »HH0,HH1 ,2 ) 

DHH0=DK2*HH0 
DHH1=DK2*HH1 
C1S2=CDK 1*SDK2 
C1C2=CDK1*CDK2 
IFICAL.LT.O. )G0 TO 20 
I F ( SAL.GT. .04 )G0 TO 20 
CNT=— 15. *CAL/SDK1 /SDK2 
CALL HANK (DK1 ,H0 ,H1 ,2 ) 

DHO= DK1 4H0 
DH1=DK1 *H1 
UKS=DK1+0K2 

CALL HANK (OKS ,SH0,SH1 ,2 ) 

DSHO=DKS*SHO 

DSH1=DKS*SH1 

Q1 1= CNT* ( C DK1 *OSH 1-C1 S2*DH0-C 1C 2*DH 1-DHH 1 + CCP*C DK2 ) 

Q12=CNT* ( CDK2=.' : DH H1-S0K2*DHH0-CCP+CDK1*DH 1 + C1 S2*DSH0-C1C2*DSH1 ) 
Q21=CNT*(SDK2*DHn-DSHl+CDK2*0H 1+CDK 1*DHH 1-CCP*C 1 C2 ) 
Q22=CNT*(C1S2*DHH0-C1C2-DHH1+CCP*CDK1~DH1-SDK2*DSH0+CDK2*DSH1) 
RETURN 

20 CONTINUE 

S 1 1 =— DHHl + CCP^CDK 2 

SI 2 =— SDK 2*DHH0+C DK2*DHH1-CCP 

S21=(0HH1-CCP*CDK2 )*CDK1 

S22 = ( SDK2-DHH0-C0K2YDHH1+CCP )*CDK1 

DKS1=DK1**2 

AL= ATAN2 < SAL ,CAL ) 

RM IN=DK1 

I F( CAL .GE .0 • ) GO TO 30 

RMIN=DK1*SAL 

DCR=— DK1- : CAL 

IFIDK2.LT .OCR ) RMI M = SORT { OKS 1+2 • *DK1*DK2$CAL+DK2*DK2 ) 

30 FNT=l + (4'+INT)/10 

INP=FNT*DK2/RMIN • 

I N P = 2 * ( INP/2) 

I F ( INP.LT.2 )1NP=2 
F I T = I N P 
I P= I NP+ 1 
DT=0K2/ F IT 
TK = . 0 

SX1=(.0,.0) 

SX2= ( .0,.0) 

SG I=-i . 


Fig. 31. Subroutine ZMM2 
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DO 90 1=1, IP 
D=SGI+3. 

IF{ I.EO.l )D=1« 

I F ( I.EO.IP )D=1. 

TKS=TK*TK 

RK= SORT (DKS1+2.*DK1#TK#CAL+TKS ) 
CALL HANK (RK,H0,H1 ,0) 

S1=SIN( DK2— TK ) 

S2=S IN ( TK ) 

SX1=SX1+S1*H0*D 
SX2=SX2+S2#HO#D 
SGI=-SGI 
TK=TK+OT 
90 CONTINUE 

SXl=SXl*DT/3. 

SX2=SX2*DT/3. 

S21=S21-SX1 

S22=S22-SX2 

SI 2= SI 2+ C OKI #5X2 

S11=S11+CDK1*SX1 

rNP=2#( INT/2I 

IP= INP+1 

FIT= INP 

JP=IP 

Tll=( ,0,.0) 

T12=( ,0,.0) 

T21= ( .0,.0) 

T22=(.0,.0) 

COT=< .0,.0) 

RS = REAL ( ZS ) 

Y11=(.0,.0) 

Y12=( .0,.0) 

Y2 1= ( *0 , *0 ) 

Y22=( .0,.0) 

B=.0 

I F ( AL.LT..05JG0 TO 210 
ALT=AL/2. 

CALT=COS(ALT) 

SALT=S IN ( ALT ) 

RCP = ( DK1+DK2 )*CALT 
RSP=( DK2-DK1 )*S AL T 
PHC=ATAN2 IRSPrRCP) 

SG I=-l • 

PH=— ALT 
DPH=AL/FIT 
DO 200 1 = 1, IP 
D=SG 1+3. 

I F < I »E0 .1 ) D= 1 « 

I F ( I .EO. IP )D=1. 

SAP=SIN(ALT+PH) 

SAM= S I N ( ALT— PH ) 

I F ( PH.LE.PHC )RMAX=0K1*SAL/SAM 

IFIPH.GT.PHC )RMAX =DK2#SAL/SAP 

DRK = RMA X/F IT 

RK=.0 

SGJ=-1. 

DT 1 1 = ( ,0,.0> 

DTI 2= ( ,0,.0) 

DT21=( .0, .0) 

DT2 2= ( .0,.0) 

DY11=(.0,.0) 

DY1 2= ( .0,.0) 


Fig. 31. (Continued) 
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DY21 = ( .0, .0) 

DY2 2 = ( .0,.0) 

RKH1 = CCP 

DO 100 J=1,JP 

C=SGJ+3. 

IF ( J.EO.l )C=1. 

I F ( J.EO.JP )C = 1. 

I F ( J.EO.l ) GO TO 94 
CALL HANK (RK,H0,H1 , 1 ) 

RKH1=RK*H1 
94 CONTINUE 

SK=RK*SAM/SAL 
TK=RK*SAP/SAL 
C1=C0S( SK ) 

C2=C0S(DK1— SK) 

S1=SIN( DK2-TK ) 

S2=S IN( TK ) 

FUN=C*RKH1 

DY11=DY1 1-FUN*C1*S1 

DY12=DY12-FUN*C1*S2 

DY21=DY21+FUN*C2*S1 

DY22=DY22+FUN*C2*S2 

SG J=-SGJ 

RK= RK+DRK 

IFIRS.LE.O. )G0 TO 100 
SS1 = S IN ( SK ) 

SS2=SIN(DK1— SK) 

DT11=DT11+FUN*SS1*S1 
DT 1 2=DT1 2+ FUN*SS1 *S2 
DT2 1 = DT 21 + FUN*SS 2*S1 
DT22=0T22+FUN*SS2*S2 
100 CONTINUE 

b=sap*drk*d 

Y11=Y11+B*DY11 
Y12=Y12+B*DY12 
Y21=Y21+B*DY21 
Y22= Y22+B*DY22 
PH=PH+DPH 
SG I=— SGI 

I F ( RS .LE.O. )G0 TO 200 
T11=T11+B*DT11 
T12=T12+B*DT12 
T21=T21+B*DT21 
T22=T22+B*DT22 
200 CONTINUE 
B=DPH/9 • 

IFfRS.GT.O. )COT = ( . 0 , 1 . ) *ZS*DPH/ ( 72. *P I*SDK 1*SDK 2*SA L) 
210 CONTINUE 

CNT=-15./SDK1/SDK2 
Q11=CNT*(CAL*S114B*Y11 )+CQT*Tll 
Q1 2=CNT* ( CAL^Sl 2+B*Yl 2 ) +COT *T12 
021=CNT*(CAL*S2l+B*Y21 )+C0T*T21 
Q22=CNT* ( CAL*S22+B*Y22 ) +COT *T22 
RETURN 
END 


Fig. 31. (Continued) 
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Fig. 32. Intersecting segments in ZMM2. 
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APPENDIX VI 
Subroutine ZMM3 


ZMM3, listed in Fig. 33, calculates the mutual impedances 
Pll, P12, P21 and P22 of two spatially-separated monopoles. In 
the input data (XI ,Y1 ) and (X2,Y2) are the coordinates of the end- 
points of segment 1. For example, XI = kx] . Similarly, (X3,Y3) and 
(X4,Y4) are the endpoints of segment 2. DK1 and DK2 denote the segment 
lengths kh-| and kh£. 

A coordinate rotation and translation is applied to move segment 
.1 onto the x axis with (XI ,Y1 ) at the origin and (X2,Y2) at (x,y) = 

(DK1 ,0) . The new coordinates of (X3,Y3) and (X4,Y4) are denoted by 
(XA,YA) and (XB,YB), respectively. Now segment 2 forms an angle a 
with the x axis, and CAL and SAL denote cos (a) and sin (a). RMIN denotes 
the shortest distance between segment 1 and segment 2 and is employed 
to determine the number of terms in the numerical integrations. 

In the DO LOOP ending at statement 50, a trigonometric table is 
calculated and stored. These sine and cosine functions represent the 
modal current (and its derivative) on segment 1. These are multiplied 
by the Simpson-rule integration coefficients before storing. 

If the cylinder has perfect conductivity, the impedance is 
given by Eqs. (68), (69) and (70). We use Eq. (24) for E x in Eq. (69). 
To obtain Ey for Eq. (70), we use a component from the second integral 
in Eq. (22). To obtain H z for the surface-impedance term aZ, we 
use Eq. (16). 

The DO LOOP ending with statement 100 uses Simpson's rule to 
integrate over segment 1 in the calculation of Ey and H z . The outer 
DO LOOP ending with statement 200 uses Simpson's rule to integrate 
J2-E-| and H] over segment 2. The impedance Z of the perfectly 
conducting cylinder is denoted by Sll , S12, etc. The term aZ is 
denoted by Til, T12, etc. 
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SUBROUTINE ZMM3 ( XI , Y1 , X2 , Y2 , X3 , Y3 , X4 , Y4 , Z S , 

2DK1,DK2, INT ,P11 ,P12,P21 ,P22 ) 

COMPLEX HHA.HHB, ZS ,HZ1 ,HZ2 ,COT ,ET 1 , E T2 , HO ,H 1 

COMPLEX P11,P12,P21,P22,S11,S12,S21,S22, T1 1, T12 , T21 , T22 

DIMENSION CC1 (21 ) » SS 1 ( 21 ) ,CC2(21 ) , SS2( 21 ) 

DATA ETA, PI/376 .7 27, 3. 14159/ 

S11=(.0,.0) 

SI 2= ( .0, .0) 

S2 1 = ( .0,.0) 

S22= ( .0, .0) 

T 1 1 = ( .0,.0) 

T12= ( .0, .0 ) 

T 2 1 = ( .0 , .0 ) 

T22= ( .0,.0) 

CBET=( X2— XI )/DKl 

SBET=( Y2-Y1 )/DKl 

XA= ( X3— XI )*CBET+(Y3-Y1 )*SBET 

XB = ( X4-X1 )*CBET+( Y4-Y1)*SBET 

YA=-< X3-X1 )*SBET+(Y3-Y1 )*CBET 

YB=-( X4-X1 )*SBET+(Y4-Y1 )*CBET 

CAL=( XB-XA)/DK2 

SAL= ( YB-YA )/DK2 

RM I N= 10000 • 

X=XA 
Y=Y A 

DX=DK2*CAL/4. 

DY=DK2*SAL/4. 

DO 40 J= 1 ,5 
YS= .0 
R=ABS( Y) 

I F ( R.GT.l.E— 15)YS=Y*Y 
XS = .0 

XAB=ABS ( X— DK 1 ) 

IF<XAB.GT.1.E-15)XS=XA8*XAB 
IF( X.LT.O. )R=SORT(X*X+YS) 

IFIX.GT.DKl )R=SORT IXS+YS ) 

I F ( R.LT .RM IN ) RM I N=R 

X=X+DX 

Y=Y+DY 

40 CONTINUE 

FNT= 1+ ( 4* I NT ) / 10 
ISS=FNT*DK1/RMIN 
ISS=2*( ISS/2 ) 

I F ( ISS.LT.2 ) ISS=2 

I F ( ISS.GT.20 ) ISS =20 

FSS= ISS 

IS0=ISS+1 

DS=DK1/FSS 

ITT = FNT *DK2/RMIN 

ITT=2*( ITT/2) 

I F ( ITT ,LT .2 ) I TT = 2 

I F ( ITT .GT .20) ITT=20 

FTT = ITT 

ITQ=ITT+1 

DT=DK2/FTT 

XP = .0 

SGN=-1. 

RS = REAL ( ZS ) 

JUM P=0 

ASAL=ABS ( SAL ) 

I F { RS .LE.0..AND. ASAL.LT..04) JUMP = 1 
IF( JUMP.EO.l )G0 TO 60 

Fig. 33. Subroutine ZMM3. 
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DO 50 1 = 1, ISO 
C= SGN+3. 

I F ( 1.60.1 )C=1 . 

IF ( I. 60. IS0)C=1. 

CC1 ( I ) =C*COS ( OKI -XP ) 

SS1( I ) =C*S IN (DK 1- XP ) 

CC2 ( I)=C*COS(XP> 

SS2( I ) =C.*S IN (XP ) 

SGN=—SGN 
XP=XP+DS 
50 CONTINUE 
60 CONTINUE 
DX=DT*CAL 
DY=DT*SAL 
X=X A 
Y= YA 
TK= .0 
SGJ=— 1 * 

CDK 1=C0S ( OKI ) 

DO 200 J = 1 * ITO 
D=SGJ+3 . 

I F ( J.EO.l ) D=1 • 

IF( J.EO .IT0)D = 1. 

CT1=D*S IN (DK2— TK) 

CT2=D*S IN ( TK ) 

XP=.0 
YS= .0 

YAB=ABS ( Y ) 

I F ( YAB.GT.l.E-15 )YS=YAB*YAB 
ET1= ( .0,.0> 

ET2 = ( .0,.0) 

HZ1=(.0,.0) 

HZ2 = ( .0, .0 ) 

RKA=SQRT (X*X+YS ) 

RKB=SQRT ( (X-DK1 ) **2+YS ) 
SPH=YAB/RKA+YAB/RKB 
I F ( SPH.LT. .04 )G0 TO 110 
I F ( JUMP. 60. 1 >G0 TO 110 
DO 100 1 = 1, ISO 
DELX=ABS ( X-XP) 

DXS= .0 

IFIDELX.GT .1 »E-15 )DXS=DElX*DELX 
RK=SORT (DXS+YS ) 

SPH=Y/RK 
C1 = CC1( I) 

S1 = SS1 ( I ) 

C2=CC2 ( I ) 

S2=SS2 ( I ) 

CALL HANK ( RK ,H0 , HI , 1 ) 
ET1=ET1-C1*SPH*H1 
ET2=ET2+C2*SPH*H1 
XP=XP+DS 

I F ( RS.LF.O. )G0 TO 100 
H Z 1 = HZ1 + S 1 *H1*SPH 
HZ2=HZ2+S2-H1^SPH 
100 CONTINUE 
110 CONTINUE 

CALL HANK ( RKA , HHA , HI ,0 ) 

CALL HANK ( RKB ,H0 ,H1 ,0 ) 
ETl=ETl*SAL*DS/3. +CAL*(CDK1*HHA-H0) 
ET2=ET2*SAL*DS/3 .+CAL* (CDK 1*H0-HHA ) 
S11=S11+CT1*ET1 


Fig. 33. (Continued) 
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S12=S12+CT2*ET1 

S21=S21+CT1*ET2 

S22 = S22+C.T2*ET2 

SGJ=-SGJ 

TK = TK + OT 

X = X+DX 

Y=Y+DY 

IFIRS.LE.O. )GO TO 200 
T11=T11+CT1#HZ1 
T12=T12+CT2*HZ1 
T21=T21+CT1*HZ2 
T22=T22+CT2*HZ2 
200 CONTINUE 

SDK1 = S IN (OKI ) 

SDK 2=S I N ( DK2 ) 

CST=-ETA*DT/ (24.*PI*SDK1$SDK2) 

COT = ( .0,1. >*DS*DT*ZS/(72.*PI*S0K1*SDK2> 

PI 1= CST*S 1 1 +COT £T 1 1 

P12=CST*S12+CQT*T12 

P2 1= CST : ^S21 +CQT *T 21 

P22=CST*S22+CQT*T22 

RETURN 

END 


Fig. 33. (Continued) 



APPENDIX VII 
Subroutine HANK 


HANK, listed in Fig. 34, generates the Hankel functions of the 
second kind H 0 (x) and H-|(x) denoted by H and HI. The input data are 
the argument X and the integer ID. Although ID is ignored in HANK, it 
is defined as follows: Let ID = 0, 1 or 2 if H 0 , H] or both are 

required. 

HANK always generates H 0 and H-] on each call. An improvement 
in computational speed could be achieved by generating just the 
required function as indicated by ID. 

In this subroutine B and B1 denote the Bessel functions and Y 
and Y1 the Neumann functions. HANK uses the polynomial approximations 
given in reference [17]. 
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SUBROUTINE HANK (X ,H .HI f 10 ) 

COMPLEX H,HI 

DATA TSP/. 63661977/ 

I F ( X.GT.3. )G0 TO 100 
XLN=TSP*AL0G(X/2. ) 

B= • 0 
B1 = .0 
Y= « 0 
Y1 = .0 
Xl= X/3 • 

X2=X1*X1 

IF( XI .LT..1 )G0 TO 60 

X4=X2*X2 

X6=X2*X4 

IF(X1.LT..3)G0 TO 55 
X8=X2*X6 
X10=X2*X8 
X12=X2*X10 

B=.21E— 34X1 2— .3 94 44E— 2*X1 0+.444479E- 1*X 8 
Y=-.24846E-34X12+.427916E-24X10-.4261214E-14X8 
Bl=. 1109 E— 44X12— .31761E -3 4X10+.443319E-24X8 
Y1=.27873E-24X12-.400976E-14X10+. 3 12 39514X8 
55 B=B-. 31638664X6+1 .26562084X4 

Y=Y+.253001174X6-. 743503844X4 
B1=B1-. 39542 89E-1 4X6+.210935734X4 
Yl= Yl-1 .31648274X6+2.16827094X4 
60 B=B-2. 24999974X2+1. 

Y=Y+. 60 559366 4X2 +. 36746691 +XLN4B 
B1=X4( B1-. 5624998 5*X2+. 5 ) 

Y1=(Y1+. 22120914 X2— .63 66198 ) /X+XLN4B1 
GO TO 200 
100 SW=SQRT (X ) 

Xl=3./X 
X2=X14X1 
X3=X14X2 
X4= XI 4X3 
X5=X1 4X4 
X6=X1*X5 

F=. 79788456-. 77E-64X1-. 55274E-2 4X2-. 951 2 E-4*X3+. 137237E-2*X4 
2-.72805E-34X5+.14476E-34X6 

T=X—. 7 85 3 98 16- .41 66397E-1 4X 1- . 3954E-44X 2+. 262 573E-2 *X3 
2-.54125E-34X4-.29333E-34X5+. 13558E-34X6 
B=F*COS ( T )/SW 
Y=F4SIN(T)/SW 

F=. 79788456+. 156F-54X1+.1659667E-14X2+. 17 105E-34X 3- . 249 51 1 E-24X4 
2+.113653E-24X5-.20033E-3*X6 

T=X-2.3561945+.124996124X1+.565E-44X2-.637879E-24X3+.74348E-3*X4 
2+.79824F-34X5— .29166E— 3*X6 
B1 = F4C0S (T )/SW 
Y1=F4SIN(T )/SW 
200 H=CMPLX ( B »— Y) 

H1 = CMPLX(B1,-Y1 ) 

RETURN 

END 


Fig. 34. Subroutine HANK. 
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APPENDIX VIII 
Subroutine CROUT 


CROUT, listed in Fig. 35, solves a system of simultaneous linear 
equations with complex coefficients. This subroutine uses the method 
of P. D. Crout described in Reference [18]. Although this subroutine 
does not use pivoting, it is efficient and accurate in the application 
discussed herein. The input data are defined as follows: 

C(I,J) Complex coefficients in the simultaneous equations 
S(I) Excitation column 

N Size of the square matrix C 

IDM Dimensions of C and S 

ISYM Zero or one for symmetric or nonsymmetric matrix 

IWR One or zero to write or suppress the solution 
112 One or two if C is original or auxiliary matrix. 

If 112 = 1, CROUT will convert the original matrix C into the 
auxiliary matrix. The auxiliary matrix is overlaid in the same 
location C, wiping out the original matrix. Similarly, the solution is 
stored in S(I) which contained the excitation column. Of course, N 
must not exceed IDM. If IWR = 1, the solution will be printed out 
with the following definitions: 

I Index number of the solution S(I) 

SNOR Normalized magnitude of S ( I ) 

SA Absolute magnitude of S(I) 

PH Phase of S(I) in degrees. 


SUBROUTINE CROUTI C ,S,N, IDM, IS YM , I WR , 1 12 ) 
COMPLEX C ( IDM , IDM) ,S U DM ) 

COMPLEX F , P »SS , T 

2 FORMAT ( IX, 1 15 ,1F10.3,1F15.7, IF 10,0 ) 

5 FORMAT ( 1 HO ) 

I F ( I12.NE.1 >G0 TO 22 
IFIN.EQ.l )S(1)=S( 1 )/C(l,l) 

I F ( N.EO.l )G0 TO 100 
I F < I SYM, NE .0 )G0 TO 8 
DO 6 1=1, N 
DO 6 J=I »N 
C< J,I)=CU»J> 

6 CONTINUE 
8 CONTINUE 

F=C( 1,1) 

DO 10 L=2,N 

Fig. 35. Subroutine CROUT. 
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10 CU,L)=C(lfU/F 
00 20 L=2 » N 
LLL=L— 1 

DO 20 I=L » N 
F=C( I,U 
DO 11 K=1,LLL 

11 F=F— C( I »K )*C (KfL) 

C ( I ,l) = F 

IF ( L «£Q, I JGO TO 20 
P=C ( L »L ) 

IF( ISYM.E0.OJGO TO 15 
F=C! L ,1 ) 

DO 12 K= 1 »LLL 

12 F=F— C (L»K)*C(K»I ) 

C ! L » I ) =F/P 

GO TO 20 
15 F=C(I,L) 

C ( L , I ) = F/P 
20 CONTINUE 
22 CONTINUE 

DO 30 1=1 ,N 
P=C!L,L) 

T=S( L) 

IF(L.E0.1)G0 TO 30 
LIL=L-1 
DO 25 K=1 » LLL 
25 T=T-C( L,K )*S(K) 

30 S<U=T/P 

DO 38 L=2,N 
I=N-L+1 
11=1+1 
T=S(I) 

DO 35 K= 1 1 »N 
35 T=T-C(I ,K)*S(K> 

38 S ( I ) =T 

IF ( IWR.LE.O) GO TO 100 
CNOR= .0 
DO 40 1=1, N 
SA=CABS(S(in 
I F ( SA.GT.CNOR )CNOR=SA 
40 CONTINUE 

I F ( CNOR «LE .0 . )CN0R=1 • 

DO 44 1=1, N 

ss=sm 

SA=CABS!S$) 

SNOR=SA /CNOR 
PH=.0 

IF( SA.GT.O. )PH=5 7.2957 8*ATAN2< AIMAGISS ),REAL(SS) ) 
WRITE! 6,2 ) I , SNOR, SA ,PH 
44 CONTINUE 
WRITE! 6,5 ) 

100 CONTINUE 
RETURN 
END 


Fig. 35. (Continued) 
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APPENDIX IX 
Subroutine VNAS 


VNAS, listed in Fig. 36, arranges the solution for a two-sided 
narrow axial slot. The excitation column CO (I) is set up to represent 
a unit voltage generator at the terminals of mode IGN. The current 
distribution induced on the cylinder is obtained by calling CROUT. 

The slot admittance Yll is then equal to the current CJ(IGN) at the 
terminals. 


SUBROUTINE VNAS ( IDM , IGN , IS YM , IWR , 1 12 ,N ,C, C J t Y 1 1 ) 
COMPLEX C( IDM, IDM ) , C J ( I DM ) , Y1 I 
DO 20 1=1, N 
20 CJ(I)=(.0,.0) 

Cj< IGN ) = ( 1 • ,0 • ) 

CALL CROUT (C ,C J ,N , I DM , I S YM, IWR, M2) 

112 = 2 

YU=CJ(IGN). 

RETURN 

END 


Fig. 36. Subroutine VNAS. 



APPENDIX X 
Subroutine VWAS 


VWAS, listed in Fig. 37, arranges the solution for a two-sided 
wide axial slot. AK denotes ka, where "a" is the aperture width. Most 
of the subroutine is concerned with generating the excitation column 
VJ(I) corresponding to a unit aperture voltage. A minor generalization 
of Eq. (38) is employed so that the aperture segments need not have 
equal lengths. The current distribution C J ( I ) induced on the cylinder 
is obtained by callinq CROUT. Finally the slot admittance Yll is cal- 
culated as in Eq. (40). 


SUBROUTINE VWAS ( I A , I B , I DM, I SYM , IWR, 1 1 , 1 2 , 1 3, 1 12 , JSA , JSB, MD , N, ND 
2,NM ,C ,CJ,D,VJ ,Y1 1) 

COMPLEX C ( I DM, I DM ) ,CJ (I DM ),VJ( I DM), Yll 

DIMENSION I A C I DM ),IB ( IDM) ,1 1 ( IDM) , 12(1 DM) ,1 3( I DM) ,MD< IDM, 5) 

2 » ND( IDM) ,D( IDM) 

AK=.0 

DO 20 K= JSA , JSB 
20 AK= AK+D (K ) 

DO 30 1 = 1, N 
30 VJ< I)=< .0, .0) 

I F < JSB.GT.JSAJGO TO 200 

K=JSA 

DK=D( K) 

V= ( l.-COS(DK) ) / ( AK*S IN ( DK ) ) 

KA= I A ( K ) 

K8= I B ( K ) 

NOK=ND(K ) 

DO 140 11=1, NDK 
I=MD(K,II) 

FI=1. 

IF ( KB «E0 • 12 ( I ) ) GO TO 136 
IFUB.EO.IKI ) ) F I=— 1 • 

GO TO 140 

136 IF( KA.E0.I3I I ) )F I=-l. 

140 VJ( I )=VJ{I)+FI*V 
GO TO 280 
200 CONTINUE 
K.A= I A( JSA ) 

KB=IB( JSA) 

LA= IA( JSA+1 ) 

LB=IB( JSA+1 ) 

IND=(LA-KB)*(LB-KB) 

I F ( IND.FQ.O)GO TO 210 
KA= I B ( JSA) 

KB=IA( JSA ) 


Fig. 37. Subroutine VWAS. 


71 



210 CONTINUE 

DO 250 K=JSA,JSB 
DK = D( K ) 

V=( 1 .-COS (DK) )/ ( A K*S I N ( OK ) ) 

NDK=ND ( K ) 

DO 240 11=1 tNDK 
I=MD(K,II ) 

F 1 = 1 . 

I F ( KB. EO. 12(1 ) )G0 TO 236 
IF(KB.EO.IKI) ) FI =- 1 • 

GO TO 240 

236 IF(KA.E0.I3(I))FI=-1. 

240 VJ{ I )=VJ( I ) +F I *V 

IFfK.EO. JSB )G0 TO 250 
L A= I A ( K+l ) 

LB= I B ( K+ 1 ) 

KA = K8 
KB=L A 

IF(LA.EO.KA)KB=LB 
250 CONTINUE 
280 CONTINUE 

DO 300 1 = 1, N 
300 C J ( I ) =VJ ( I ) 

CALL CROUT (C ,C J ,N , IDM , IS YM , IWR , I 12 ) 
112 = 2 

VI 1= ( .0,.0) 

DO 400 1=1, N 

400 Y11 = Y11 + VJ(I )*CJ( I ) 

RETURN 

END 


Fig. 37. (Continued) 
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APPENDIX XI 
Subroutine VMLS 


VMLS, listed in Fig. 38, arranges the solution for a parallel 
magnetic line source near a cylinder. Most of the subroutine is con- 
cerned with generating the excitation column for a one-volt magnetic 
line source. The integrations. in Eq. (31) are performed in CMLS. The 
current distribution C J ( I ) induced on the cylinder is obtained by 
calling CROUT . Finally the admittance Vll of the line source is cal- 
culated. 


SUBROUTINE VMLS< IA , IB , IDM , I NT , I SYM , I WR , 1 1 , 1 2 , I 3 , 1 12, MD » N, ND, NM , 
2C,CJ ,D,PSI ,VJ,X ,Y ,XS,YS,Y11 ,ZS) 

COMPLEX C( IDM, IDM) ,CJ ( IDM) ,VJ(IDM) ,Y11,P1,P2, 01,02, ZS 
DIMENSION IA(IDM) , I B ( IDM ) , 1 1 ( I DM ) , I 2 ( I DM ) , I 3 ( I DM) 

DIMENSION MD( IDM ,5 ),ND < IDM ) ,X ( IDM) ,YUDM) ,D( IDM) 

DATA ETA, TP/376. 727, 6. 28318/ 

00 100 1=1, N 
VJ( I)=(.0,.0) 

100 CJ( I)=( .0,.0) 

DO 240 K = 1 » NM 
KA= IA ( K ) 

KB= I B ( K ) 

CALL CMLS (PS I ,X< KA ) , Y ( KA ) , X ( KB ) , Y ( KB ) ,XS , YS , D ( K ) , I NT , PI , P2 , Q1 , 02 ) 
s 01=ZS*01 
Q2= ZS*Q2 
NDK=ND(K ) 

DO 240 11=1, NDK 

1 = MD(K, I I ) 

F 1= 1 • 

IFIKB.EO. 12(D) GO TO 236 
IFIKB.EO.IKI )) F I =— 1 • 

CJ( I ) =CJ ( I )+FI*Pl 
VJ( I ) =VJ { I )+FI*( Pl+Ql) 

GO TO 240 

236 IFIKA.EO. 13(D) FI=-1. 

CJ( I ) =CJ ( I )+FI*P2 
VJ( I)=VJ(I )+FPM P2+02) 

240 CONTINUE 

CALL CROUT (C »CJ ,N,IDM, ISYI-1, 1 V'R ,112) 

112=2 

Y11=CMPLX(TP/ (4. *ETA ) ,0. ) 

DO 300 1=1, N 
300 Y11=Y11+CJ ( I )*VJ (I ) 

RETURN 

END 


Fig. 38. Subroutine VMLS. 
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APPENDIX XII 
Subroutine CMLS 


CMLS, listed in Fig. 39, evaluates the coupling between a 
one-volt magnetic line source and an electric strip monopole. The 
line source has coordinates XS and YS. The monopole (or segment) has 
length DK and endpoint coordinates (X1,Y1) and (X2,Y2). PI is the 
coupling to the mode with terminals at point 1, and P2 applies to 
terminals at point 2. The excitation voltage V m in Eq. (31) is 
obtained by integrating over both segments of test-dipole m. CMLS 
integrates over one segment to obtain PI and P2, and VMLS adds the 
appropriate quantities (PI or P2) from two segments to obtain V m . 

D1 and D2 are the distances between the line source and (XI ,Y1 ) 
and (X2,Y2) respectively. RMIN is the shortest distance between the 
line source and the segment. 

A coordinate rotation and translation moves the segment onto 
the x axis with (XI ,Y1 ) at the origin and (X2,Y2) at x = DK. The 
new coordinates of the line source are (XA,YA). The integration in 
Eq. (31) is evaluated with the trapezoidal rule in the DO LOOP ending 
at statement 100. Eq. (41) is employed for the electric field of the 
magnetic line source. 

The admittance of a one- volt magnetic line source near a 
cylinder is defined as the magnetic field intensity H z at the line 
source. The field H z is the sum of the free-space field of the line 
source, the field of the electric current on the cylinder, and the 
field of the magnetic current on the cyliner. 

By reciprocity, PI and P2 are useful not only in generating the 
excitation column but also in calculating H^. To permit VMLS to 
calculate the line-source admittance, CMLS also generates Q1 and Q2 
which denote the field from the magnetic current on one segment 
of the cylinder. Equations (6) and (42) are employed here, and the 
integrations for Q1 and Q2 are performed in the DO LOOP ending with 
statement 100. 


74 



SUBROUTINE CML S ( P S I ,X 1 , Y 1 ,X2 , Y2 ,XS , YS , OK , I NT t P 1 , P2 , Q 1 , 02 ) 
COMPLEX CST,H0,H1,P1,P2,Q1,Q2 
DATA ETA/376.727/ 

DKH=DK/ 25 • 

D I = SQRT ( (XS-X1 )**2+( YS-Y1 )**?.) 

D2= SORT ( (XS-X2 >**2 + < YS-Y2 )**2 ) 

P2=(.0,.0) 

P1=CMPLX( .5*PS 1/360., 0. ) 

01 =( . 0 ,. 0 ) 

02 =( . 0 , . 0 ) 

IFfDl.LT. DKH )G0 TO 200 
P1=(.0,.0) 

P2=CMPLX ( .5*PS 1/360. ,0. ) 

I F ( D2 .LT .DKH ) GO TO 200 
SDK= S I N ( DK ) 

P1=(.0,.0) 

P2 = ( .0,.0) 

CBET= ( X2-X1 ) / DK 

SBET = ( Y2-Y1 ) /DK 

XA= (XS-X1 )*CBET+ ( YS-YI )*SBET 

YA=— ( XS— XI )*SBET+ ( YS-Y1 )*CBET 

X=X A 

Y= YA 

YSQ=Y**2 
RM IN=ABS f Y ) 

I F ( X.LT .0. )RMIN=SORT(X*X+YSQ) 

IFIX.GT.DK )RMIN=SORT( ( X-DK ) **2 +YSQ ) 

FNT = 1+ ( 4* I NT ) / 10 
I SS= FNT *DK/ RM IN 
I F ( ISS.LT.2)ISS=2 
F IT= ISS 
DS=DK/F IT 
XP=DS/2. 

DO 100 1=1, ISS 
DELX=X— XP 

RK=SORT (DELX*~2+YSQ) 

SPH=Y/RK 
S1=SIN( DK-XP) 

S2 = S IN ( XP ) 

CALL HA NK ( RK , HO , HI , 2 ) 

P1=P1+S1 *H1*SPH 
P2~P2+S2*H1*SPH 
Q1=01+S1*H0 
Q2=Q2+S2*HO 
XP=XP+DS 
100 CONTINUE 

CST= ( .0,1. )*DS/ (4 ,*SDK) 

P1=CST*P1 

P2=CST*P2 

CNT=-DS/ (4 .*ETA* SDK ) 

01=CNT *01 
02=CNT *02 
200 CONTINUE 
RETURN 
END 


Fig. 39. Subroutine CMLS. 
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APPENDIX XIII 
Subroutine VFF 


VFF is listed in Fig. 40. In antenna problems (options 1, 2 
and 3) and bistatic scattering (option 4 with INC = 0), VFF calculates 
the magnetic field HZT in the far-zone and the GAIN or the echo width 
EWL. In these cases, the current distribution CJ ( I ) on the cylinder 
is already known. VFF calls CFF to obtain the far-field from each 
strip on the cylinder. HJJ(I) and HMM(I) denote the fields from the 
electric and magnetic mode currents, respectively, of mode I. HZS is 
the field generated by the currents on the cylinder, HZM is the field 
of the magnetic line source in option 3, and HZT is the sum of these. 

In backscattering problems (option 4 with INC = 1 ) , VFF 
generates the excitation column which is related to HJJ(I) by the 
reciprocity theorem. Then VFF calls CROUT to obtain the current 
distribution C J ( I ) and calculates the extinction cross-section ECS, 
the backscattered field and the echo width. 
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SUBROUTINE VFF ( I A , IB , I NC , I DM , I SYM , IWR , 1 1 , I 2 , I 3 , 1 1 2 , LOP, MD , N, ND , NM , 
2C,CJ,D,EWL ,G,GAIM ,HJJ ,HMM ,HZS ,HZT,PH,ECS, VJ,X ,Y,XS, YS,ZS) 

COMPLEX COT , DOT , HJ 1 , HJ 2 ,HM 1 « HM2 ,HZM,HZS,HZT ,ZS 
COMPLEX C J ( I DM ) ,H J J ( I DM ) , HMM ( I DM ) ,C ( I DM, I DM ) , VJ ( I DM ) 

DIMENSION I A ( I DM ), IB ( I DM) , I 1 ( I DM ) , 12 ( I DM) , I 3 ( I DM) , MD ( I DM, 5) 
DIMENSION ND(IDM) ,X ( I DM ) , Y ( I DM ) , D ( I DM ) 

DATA ETA, TP/376. 727,6.28318/ 

COT=l .414214*ETA*CMPLX(1. ,-l. ) 

I F ( IS YM .NE .0 ) DOT =CQT*CONJG!ZS)/ZS 
EC S= .0 

PHR=.0174533*PH 
C PH=COS ( PHR ) 

SPH = S IN ( PHR ) 

DO 232 1=1, N 

H J J ( I ) = ( .0 , .0 ) 

232 HMM ( I )=( .0, .0) 

DO 250 K=1 »NM 
KA= IA( K ) 

K8= IB ( K ) 

CALL CFF(X(KA),Y(KA),X(KB ),Y(KR),D(K) 

2 , CPH , SPH , ZS ,H J 1 , HJ2 ,HM 1 ,HM2 ) 

NDK= ND ( K ) 

DO 240 11=1, NDK 
I = MD ( K , I I ) 

F 1= 1 . 

IF(KB.E0.I2( I ) ) GO TO 236 
I F ( KB.EO.il (I) )F I=-l . 

H J J ( I)=HJJ(I)+FI*HJ1 
HMM ( I ) =HMM ( I ) +FI *HM1 
GO TO 240 

236 1 F ( KA.E0.I3I I ) )F I=-l. 

HJJ( I)=HJJ(I)+FI*HJ2 
HMM ( I ) =HMM ( I ) +F I *HM2 
240 CONTINUE 
250 CONTINUE 

I F ( INC .L E .0 ) GO TO 270 
DO 260 1=1, N 
CJ< I )=COT*HJJCI } 

VJ ( I ) =C J ( I > 

260 I F ( ISYM.NE.O )VJ (I )=VJ(I ) -DOT #HMM { I ) 

CALL CROUT (C,CJ,N,IDM, ISYM , I WR, 1 12 ) 

112=2 

DO 265 1=1 , N 

265 ECS=ECS+REAL (VJ (I )*CONJG (CJ( I ) ) ) 

ECS=ECS/ETA 
270 HZS= ( .0,.0) 

DO 360 1=1, N 

360 HZS=HZS+C J ( I )* ( HJ J ( I ) +HMM ( I ) ) 

HAB = CAB S ( HZS ) 

IF ( LOP .E0.4 )EWL=T P*HA8*HAB 
HZT = HZS 

IF(L0P.ME.3)G0 TO 400 

PS I = XS YCPH+YS *SP H 

HZM=CMPLX (COS ( PST ),SIN(PSI)) 

HZM=-( l.,l . )*HZM/(2.*1 .4142 14*ETA ) 

HZT= HZS+HZM 
400 HAB=CABS(HZT ) 

IF(L0P.LT.4)GAIN=TP*ETA*HAB*HAB/G 

RETURN 

END 

Fig. 40. Subroutine VFF. 
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APPENDIX XIV 
Subroutine CFF 


CFF, listed in Fig. 41, calculates the far-zone field of a strip 
monopole. The endpoints of the monopole are at (XA,YA) and (XB,YB), and 
DK is the length. CPH and SPH denote cos<j> and sin<j>. HJ1 and HJ2 are 
the far-field contributions from the electric mode currents with 
terminals at (XA,YA) and (XB,YB) respectively. Similarly HM1 and HM2 
are the fields of the magnetic mode currents. This subroutine uses 
Eqs. (44) through (57). 


SUB ROUT INE CFF ( X A , Y A , XB , VB , DK ,CPH , SPH , ZS , H J L ,H J 2 , HI-11 , HM2) 
COMPLEX EJA,EJB ,C ST , ZS , H J 1 , H.I2 , HM 1 , HM2 
DATA ETA, PI/376. 727,3.14159/ 

CA= ( XB-XA )/DK 
CB=< YB-YA )/DK 

G=CA*CPH+CB*SPH 

P=CB-CPH-CA-SPH 

GK=P**2 

A = X A*C P H+ YA^S PH 
B=XB*CPH+YB*SPH 
E JA = CMPLX { COS ( A ) t S IN ( A ) ) 

EJB=CMPLX (COS (B ), SIN(B) ) 

SDK = S IN (DK ) 

CDK=COS ( DK ) 

IF( GK.LT..001 )GO TO 250 

CST=CMPLX (1. ,1. )/ (4.*PI*SDK*1.414214*GK) 

HM1=CST* ( EJA*CMPLX (CDK ,G*SDK )-EJB) 

HM2=CST* ( E JB*CMPL X ( CDK , — G*S DK } - E J A ) 

GO TO 300 

250 CST=CMPIX (-1. ,1 .) / ( 8 . *P I *1 . 4 1 42 14*SDK ) 

IF( G.LT *0. )G0 TO 280 
HM1=CST*(DK*EJB~SDK*EJA ) 

HM2=CST * ( SDK*E JB-DK^EJ A ) 

GO TO 300 

280 HM1=CST*(SDK*EJA-DK*EJB) 

HM2=CST* ( DKSEJA-SDK4EJB ) 

300 HJ1=P*HM1 - " - - ' • 

HJ 2= P*HM2 

HM1=-ZS*HM1/ETA 

HM2=— ZS*HM2/ETA 

RETURN 

END 


Fig. 41. Subroutine CFF. 
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APPENDIX XV 
Subroutine C SURF 


CSURF, listed in Fig. 42, calculates the surface impedance ZS 
of a plane conducting slab for normal incidence. FMC denotes the 
frequency in MHz. The conductivity of the slab material in megamhos 
per meter is CMM. The slab thickness is t and TK denotes kt. 


SUBROUTINE CSURF (CMM, FMC ,TK,ZS ) 

COMPLEX ETA,R,ZS» ETBT 

DATA E»FTA0»TP»U/8. 8 54 33E-12»376*7 27 ,6*28318,12.5 664 E— 7/ 
ALPH=SORT (TR*FMC*U*CMM/2. ) * 1 . E 6 
SOT=SORT(TP*FMC*U/ (2.*CMM) ) 

£T A=CMPL X (SOT , S OT I 

TAT = 2 .*TK*SQRT ( C MM/ ( 2 . *E*T P*FMC ) ) 

ZS=ETA 

I F ( TAT.GT.60. )GO TO 100 
ETAT=EXP (-TAT ) 

ETBT=CMPLX (COS (T AT ) ,-S IN (TAT ) ) 

R=ETAT*ETBT*(ETAO-ETA )/ (ETAO+ETA) 

ZS=ETA*(1.+R)/(1 .-R) 

100 continue 

RETURN 

END 


Fig. 42. Subroutine CSURF. 
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